ThreeB 1.1

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

- File List
- File Members

sampled_multispot.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 SPOT_WITH_BACKGROUND_H
    00021 #define SPOT_WITH_BACKGROUND_H
    00022
    00023 #include <vector>
    00024 #include <cvd/image_ref.h>
    00025 #include <tr1/tuple>
    00026 #include <TooN/TooN.h>
    00027
    00028 #include "drift.h"
    00029
    00030 typedef char State;
    00031
    00032 namespace SampledMultispot
    00033 {
    00034
    00035 using namespace std;
    00036 using namespace CVD;
    00037 using namespace TooN;
    00038 using namespace std::tr1;
    00039
    00040
    00041
    00042 //The changes to SpotWithBackground to operate with/without drift and masking are
    00043 //irritatingly pervasive. I can't think of any other way of doing it.
    00044 #define SWBG_NAME SpotWithBackground
    00045 #include "spot_with_background.hh"
    00046
    00047
    00048 #define SWBG_NAME SpotWithBackgroundMasked
    00049 #define SWBG_HAVE_MASK
    00050 #include "spot_with_background.hh"
    00051
    00052 //#define SWBG_NAME SpotWithBackgroundDrift
    00053 //#define SWBG_HAVE_DRIFT
    00054 //#include "spot_with_background.hh"
    00055
    00056 //#undef SWBG_NAME
    00057 //#define SWBG_NAME SpotWithBackgroundDriftMasked
    00058 //#define SWBG_HAVE_DRIFT
    00059 //#define SWBG_HAVE_MASK
    00060 //#include "spot_with_background.hh"
    00061
    00062
    00063
    00064
    00065 inline double intensity(double i)
    00066 {
    00067     return i;
    00068 }
    00069
    00070 inline double intensity(const pair<double, Vector<4> >& i)
    00071 {
    00072     return i.first;
    00073 }
    00074
    00075
    00076 //Add and remove a spot over the entire region
    00077 template<class T>
    00078 void remove_spot(vector<vector<double> >& current_sample_intensities, const vector<T>& spot_intensities, const vector<State>& spot_sample)
    00079 {
    00080     for(unsigned int frame=0; frame < current_sample_intensities.size(); frame++)
    00081         if(spot_sample[frame] == 0) //Spot is on, so remove it
    00082             for(unsigned int p=0; p < spot_intensities.size(); p++)
    00083                 current_sample_intensities[frame][p] -= intensity(spot_intensities[p]);
    00084 }
    00085
    00086 template<class T>
    00087 void add_spot(vector<vector<double> >& current_sample_intensities, const vector<T>& spot_intensities, const vector<State>& spot_sample)
    00088 {
    00089     for(unsigned int frame=0; frame < current_sample_intensities.size(); frame++)
    00090         if(spot_sample[frame] == 0) //Spot is on, so add it
    00091             for(unsigned int p=0; p < spot_intensities.size(); p++)
    00092                 current_sample_intensities[frame][p] += intensity(spot_intensities[p]);
    00093 }
    00094
    00095
    00096 //Add and remove a spot only over a mask.
    00097 template<class T>
    00098 void remove_spot(vector<vector<double> >& current_sample_intensities, const vector<T>& spot_intensities, const vector<State>& spot_sample, const vector<int>& mask)
    00099 {
    00100     for(unsigned int frame=0; frame < current_sample_intensities.size(); frame++)
    00101         if(spot_sample[frame] == 0) //Spot is on, so remove it
    00102             for(unsigned int p=0; p < mask.size(); p++)
    00103                 current_sample_intensities[frame][mask[p]] -= intensity(spot_intensities[mask[p]]);
    00104 }
    00105
    00106 template<class T>
    00107 void add_spot(vector<vector<double> >& current_sample_intensities, const vector<T>& spot_intensities, const vector<State>& spot_sample, const vector<int>& mask)
    00108 {
    00109     for(unsigned int frame=0; frame < current_sample_intensities.size(); frame++)
    00110         if(spot_sample[frame] == 0) //Spot is on, so add it
    00111             for(unsigned int p=0; p < mask.size(); p++)
    00112                 current_sample_intensities[frame][mask[p]] += intensity(spot_intensities[mask[p]]);
    00113 }
    00114
    00115
    00116 //Add and remove a drifty spot only over a mask.
    00117 template<class T>
    00118 void remove_spot(vector<vector<double> >& current_sample_intensities, const vector<vector<T> > & spot_intensities, const vector<State>& spot_sample, const vector<int>& mask)
    00119 {
    00120     const int steps = spot_intensities.size();
    00121     const int frames = current_sample_intensities.size();
    00122
    00123     for(int frame=0; frame < frames; frame++)
    00124     {
    00125         int s = frame * steps / frames;
    00126
    00127         if(spot_sample[frame] == 0) //Spot is on, so remove it
    00128             for(unsigned int p=0; p < mask.size(); p++)
    00129                 current_sample_intensities[frame][mask[p]] -= intensity(spot_intensities[s][mask[p]]);
    00130     }
    00131 }
    00132
    00133 template<class T>
    00134 void add_spot(vector<vector<double> >& current_sample_intensities, const vector<vector<T> >& spot_intensities, const vector<State>& spot_sample, const vector<int>& mask)
    00135 {
    00136     const int steps = spot_intensities.size();
    00137     const int frames = current_sample_intensities.size();
    00138
    00139     for(int frame=0; frame < frames; frame++)
    00140     {
    00141         int s = frame * steps / frames;
    00142
    00143         if(spot_sample[frame] == 0) //Spot is on, so add it
    00144             for(unsigned int p=0; p < mask.size(); p++)
    00145                 current_sample_intensities[frame][mask[p]] += intensity(spot_intensities[s][mask[p]]);
    00146     }
    00147 }
    00148
    00149
    00150
    00151 //Compute the spot intensity for a given spot at each pixel
    00152 inline vector<double> compute_spot_intensity(const vector<ImageRef>& pixels, const Vector<4>& params)
    00153 {
    00154     vector<double> intensities(pixels.size());
    00155
    00156     for(unsigned int i=0; i < pixels.size(); i++)
    00157         intensities[i] = spot_shape(vec(pixels[i]), params);
    00158
    00159     return intensities;
    00160 }
    00161
    00162 //Compute the spot intensity derivatives for a given spot at each pixel
    00163 inline vector<pair<double, Vector<4> > > compute_spot_intensity_derivatives(const vector<ImageRef>& pixels, const Vector<4>& params)
    00164 {
    00165     vector<pair<double, Vector<4> > > derivatives(pixels.size());
    00166
    00167     for(unsigned int i=0; i < pixels.size(); i++)
    00168         derivatives[i] = spot_shape_diff_position(vec(pixels[i]), params);
    00169     return derivatives;
    00170 }
    00171
    00172 inline vector<tuple<double, Vector<4>, Matrix<4> > > compute_spot_intensity_hessian(const vector<ImageRef>& pixels, const Vector<4>& params)
    00173 {
    00174     vector<tuple<double, Vector<4>, Matrix<4> > > hessian(pixels.size());
    00175
    00176     for(unsigned int i=0; i < pixels.size(); i++)
    00177         hessian[i] = spot_shape_hess_position(vec(pixels[i]), params);
    00178     return hessian;
    00179 }
    00180
    00181
    00182 /**
    00183 Create a sequence of integers. These can be used as observations
    00184 in an observation class by forward_algorithm() and etc.
    00185 @param n Length of sequence
    00186 @ingroup gUtility
    00187 */
    00188 inline vector<int> sequence(int n)
    00189 {
    00190     vector<int> v;
    00191     for(int i=0; i < n; i++)
    00192         v.push_back(i);
    00193     return v;
    00194 }
    00195
    00196 /*struct RndGrand48
    00197 {
    00198     double operator()()
    00199     {
    00200         return drand48();
    00201     }
    00202 };*/
    00203
    00204 ///Draw samples from the spot states given the spots positions and some data.
    00205 ///Variable naming matches that in FitSpots.
    00206 ///@ingroup gStorm
    00207 class GibbsSampler
    00208 {
    00209     const vector<vector<double> >& pixel_intensities;
    00210     const vector<vector<double> >& spot_intensities;
    00211     const vector<Vector<4> > spots;
    00212     const Matrix<3> A;
    00213     const Vector<3> pi;
    00214     const double base_variance;
    00215     double variance;
    00216
    00217     const int sample_iterations;
    00218     const int num_frames, num_pixels;
    00219     const vector<int> O;
    00220
    00221     vector<vector<State> > current_sample;
    00222     vector<vector<double> > current_sample_intensities;
    00223
    00224     public:
    00225
    00226     GibbsSampler(const vector<vector<double> >& pixel_intensities_,
    00227                  const vector<vector<double> >& spot_intensities_,
    00228                  const vector<Vector<4> >& spots_,
    00229                  const Matrix<3> A_,
    00230                  const Vector<3> pi_,
    00231                  double variance_,
    00232                  int sample_iterations_)
    00233     :pixel_intensities(pixel_intensities_), //pixel_intensities: [frame][pixels]
    00234      spot_intensities(spot_intensities_),   //spot_intensities: [spot][pixel]
    00235      spots(spots_),
    00236      A(A_),
    00237      pi(pi_),
    00238      base_variance(variance_),
    00239      variance(variance_),
    00240      sample_iterations(sample_iterations_),
    00241      num_frames(pixel_intensities.size()),
    00242      num_pixels(pixel_intensities[0].size()),
    00243         //Observations vector. As usual for this application, the observations are just
    00244         //numbered integers which refer to data held elsewhere.
    00245      O(sequence(num_frames)),
    00246         //Start all spots OFF, so the intensity is 0. OFF is 1 or 2, not 0!!!
    00247         //sample_list: [sample][spot][frame]: list of samples drawn using Gibbs sampling
    00248      current_sample(spots.size(), vector<State>(num_frames, 2)), //current sample [spot][frame]
    00249         //pixel intensities assosciated with the current sample [frame][pixel]
    00250      current_sample_intensities(num_frames, vector<double>(num_pixels))
    00251     {
    00252         //Check a bunch of stuff
    00253         assert_same_size(pixel_intensities);
    00254         assert_same_size(spot_intensities);
    00255
    00256     }
    00257
    00258     ///Update the noide variance. Used for adding thermal noise.
    00259     ///@param v noise variance.
    00260     void set_variance(double v)
    00261     {
    00262         variance = v;
    00263     }
    00264
    00265
    00266     ///Reset the gibbs sampler oro the initial state (all spots off)
    00267     void reset()
    00268     {
    00269         vector<State> off(num_frames, 2);
    00270         fill(current_sample.begin(), current_sample.end(), off);
    00271
    00272         vector<double> black(num_pixels);
    00273         fill(current_sample_intensities.begin(), current_sample_intensities.end(), black);
    00274         variance = base_variance;
    00275     }
    00276
    00277     ///Get the next sample
    00278     ///@param rng Random number generator
    00279     template<class T> void next(T& rng)
    00280     {
    00281         for(int j=0; j < sample_iterations; j++)
    00282             for(int k=0; k < (int) spots.size(); k++)
    00283             {
    00284                 //Subtract off the spot we're interested in.
    00285                 remove_spot(current_sample_intensities, spot_intensities[k], current_sample[k]);
    00286
    00287                 //Now current_sample_intensities is the image value for every spot in every frame,
    00288                 //except the current spot, which is always set to off. This allows us to add it in
    00289                 //easily.
    00290                 SpotWithBackground B(current_sample_intensities, spot_intensities[k], pixel_intensities, variance);
    00291                 vector<array<double, 3> > delta = forward_algorithm_delta(A, pi, B, O);
    00292                 current_sample[k] = backward_sampling<3,State, T>(A, delta, rng);
    00293
    00294                 //Put the newly sampled spot in
    00295                 add_spot(current_sample_intensities, spot_intensities[k], current_sample[k]);
    00296             }
    00297     }
    00298     /*void next()
    00299     {
    00300         RngDrand48 rng;
    00301         next(rng);
    00302     }*/
    00303
    00304     ///Retrieve the current sample
    00305     const vector<vector<State> >& sample() const
    00306     {
    00307         return current_sample;
    00308     }
    00309     ///Retrieve the intensities for the current sample
    00310     const vector<vector<double> >& sample_intensities() const
    00311     {
    00312         return current_sample_intensities;
    00313     }
    00314
    00315 };
    00316
    00317 ///Gibbs sampling class which masks spots to reduce computation.
    00318 ///
    00319 ///This draws samples from, the spot states given the spots positions and some data. It is
    00320 ///very similar to GibbsSampler, except that it only computes probabilities in a mask around each spot
    00321 ///to save on computation. Variable naming matches that in FitSpots.
    00322 ///@ingroup gStorm
    00323 class GibbsSampler2
    00324 {
    00325     const vector<vector<double> >& pixel_intensities;
    00326     const vector<vector<double> >& spot_intensities;
    00327     const vector<Vector<4> > spots;
    00328     const std::vector<std::vector<int> >& spot_pixels;
    00329     const Matrix<3> A;
    00330     const Vector<3> pi;
    00331     const double base_variance;
    00332     double variance;
    00333
    00334     const int sample_iterations;
    00335     const int num_frames, num_pixels;
    00336     const vector<int> O;
    00337
    00338     vector<vector<State> > current_sample;
    00339     vector<vector<double> > current_sample_intensities;
    00340
    00341     vector<double> cutout_spot_intensities;
    00342     vector<vector<double> > cutout_pixel_intensities;
    00343     vector<vector<double> > cutout_current_sample_intensities;
    00344
    00345     public:
    00346
    00347     GibbsSampler2(const vector<vector<double> >& pixel_intensities_,
    00348                  const vector<vector<double> >& spot_intensities_,
    00349                  const vector<Vector<4> >& spots_,
    00350                  const vector<vector<int> >& spot_pixels_,
    00351                  const Matrix<3> A_,
    00352                  const Vector<3> pi_,
    00353                  double variance_,
    00354                  int sample_iterations_)
    00355     :pixel_intensities(pixel_intensities_), //pixel_intensities: [frame][pixels]
    00356      spot_intensities(spot_intensities_),   //spot_intensities: [spot][pixel]
    00357      spots(spots_),
    00358      spot_pixels(spot_pixels_),
    00359      A(A_),
    00360      pi(pi_),
    00361      base_variance(variance_),
    00362      variance(variance_),
    00363      sample_iterations(sample_iterations_),
    00364      num_frames(pixel_intensities.size()),
    00365      num_pixels(pixel_intensities[0].size()),
    00366         //Observations vector. As usual for this application, the observations are just
    00367         //numbered integers which refer to data held elsewhere.
    00368      O(sequence(num_frames)),
    00369         //Start all spots OFF, so the intensity is 0. OFF is 1 or 2, not 0!!!
    00370         //sample_list: [sample][spot][frame]: list of samples drawn using Gibbs sampling
    00371      current_sample(spots.size(), vector<State>(num_frames, 2)), //current sample [spot][frame]
    00372         //pixel intensities assosciated with the current sample [frame][pixel]
    00373      current_sample_intensities(num_frames, vector<double>(num_pixels)),
    00374      cutout_pixel_intensities(num_frames),
    00375      cutout_current_sample_intensities(num_frames)
    00376     {
    00377         //Check a bunch of stuff
    00378         assert_same_size(pixel_intensities);
    00379         assert_same_size(spot_intensities);
    00380     }
    00381
    00382     ///Update the noide variance. Used for adding thermal noise.
    00383     ///@param v noise variance.
    00384     void set_variance(double v)
    00385     {
    00386         variance = v;
    00387     }
    00388
    00389
    00390     ///Reset the gibbs sampler oro the initial state (all spots off)
    00391     void reset()
    00392     {
    00393         vector<State> off(num_frames, 2);
    00394         fill(current_sample.begin(), current_sample.end(), off);
    00395
    00396         vector<double> black(num_pixels);
    00397         fill(current_sample_intensities.begin(), current_sample_intensities.end(), black);
    00398         variance = base_variance;
    00399     }
    00400
    00401     ///Get the next sample
    00402     ///@param rng Random number generator
    00403     template<class T> void next(T& rng)
    00404     {
    00405
    00406 //double remove=0;
    00407 //double cut=0;
    00408 //double swb=0;
    00409 //double ff_masked=0;
    00410 //double bs=0;
    00411 //double add=0;
    00412 //cvd_timer t;
    00413     std::vector<array<double, 3> > delta3;
    00414         for(int j=0; j < sample_iterations; j++)
    00415             for(int k=0; k < (int) spots.size(); k++)
    00416             {
    00417 //t.reset();
    00418                 //Subtract off the spot we're interested in.
    00419                 remove_spot(current_sample_intensities, spot_intensities[k], current_sample[k], spot_pixels[k]);
    00420 //remove+=t.reset();
    00421 /*
    00422                 //Cut out
    00423                 //spot
    00424                 cutout_spot_intensities.resize(spot_pixels[k].size());
    00425                 for(unsigned int i=0; i < spot_pixels[k].size(); i++)
    00426                     cutout_spot_intensities[i] = spot_intensities[k][spot_pixels[k][i]];
    00427
    00428                 //others
    00429                 for(int f=0; f < num_frames; f++)
    00430                 {
    00431                     cutout_current_sample_intensities[f].resize(spot_pixels[k].size());
    00432                     cutout_pixel_intensities[f].resize(spot_pixels[k].size());
    00433                     for(unsigned int i=0; i < spot_pixels[k].size();i++)
    00434                     {
    00435                         cutout_current_sample_intensities[f][i] = current_sample_intensities[f][spot_pixels[k][i]];
    00436                         cutout_pixel_intensities[f][i] = pixel_intensities[f][spot_pixels[k][i]];
    00437                     }
    00438                 }*/
    00439 //cut += t.reset();
    00440                 //Now current_sample_intensities is the image value for every spot in every frame,
    00441                 //except the current spot, which is always set to off. This allows us to add it in
    00442                 //easily.
    00443
    00444 //              SpotWithBackground B(current_sample_intensities, spot_intensities[k], pixel_intensities, variance);
    00445 //              vector<array<double, 3> > delta = forward_algorithm_delta(A, pi, B, O);
    00446
    00447 //ff+=t.reset();
    00448
    00449 //              SpotWithBackground B2(cutout_current_sample_intensities, cutout_spot_intensities, cutout_pixel_intensities, variance);
    00450 //              std::vector<array<double, 3> > delta2 = forward_algorithm_delta(A, pi, B2, O);
    00451 //ff_cut+=t.reset();
    00452
    00453                 SpotWithBackgroundMasked B3(current_sample_intensities, spot_intensities[k], pixel_intensities, variance, spot_pixels[k]);
    00454 //swb += t.reset();
    00455                 forward_algorithm_delta2<3>(A, pi, B3, O, delta3);
    00456 //f_masked+=t.reset();
    00457                 /*for(unsigned int i=0; i < delta.size(); i++)
    00458                 {
    00459                     cout.precision(20);
    00460                     cout.setf(cout.scientific);
    00461                     std::cout << delta[i][0] << " " << delta[i][1] << " " <<delta[i][2] << std::endl;
    00462                     std::cout << delta2[i][0] << " " << delta2[i][1] << " " <<delta2[i][2] << std::endl;
    00463                     cout << endl;
    00464                 }
    00465                 std::exit(1);
    00466
    00467 */
    00468
    00469                 current_sample[k] = backward_sampling<3,State, T>(A, delta3, rng);
    00470 //bs += t.reset();
    00471                 //Put the newly sampled spot in
    00472                 add_spot(current_sample_intensities, spot_intensities[k], current_sample[k], spot_pixels[k]);
    00473 //add += t.reset();
    00474             }
    00475 //      cout << "remove=" <<remove << " cut=" << cut << " swb=" << swb<< " ff_mask=" << ff_masked << " bs=" <<bs << " add="<<add << endl;
    00476     }
    00477 /*  void next()
    00478     {
    00479         RngDrand48 rng;
    00480         next(rng);
    00481     }
    00482 */
    00483     ///Retrieve the current sample
    00484     const vector<vector<State> >& sample() const
    00485     {
    00486         return current_sample;
    00487     }
    00488     ///Retrieve the intensities for the current sample
    00489     const vector<vector<double> >& sample_intensities() const
    00490     {
    00491         return current_sample_intensities;
    00492     }
    00493
    00494 };
    00495
    00496 #if 0
    00497 ///Gibbs sampling class
    00498 class GibbsSampler3
    00499 {
    00500     const vector<vector<double> >& pixel_intensities;
    00501     const vector<vector<vector<double> > >& spot_intensities;
    00502     const vector<Vector<4> > spots;
    00503     const std::vector<std::vector<int> >& spot_pixels;
    00504     const Matrix<3> A;
    00505     const Vector<3> pi;
    00506     const double base_variance;
    00507     double variance;
    00508
    00509     const int sample_iterations;
    00510     const int num_frames, num_pixels;
    00511     const vector<int> O;
    00512
    00513     vector<vector<State> > current_sample;
    00514     vector<vector<double> > current_sample_intensities;
    00515
    00516     public:
    00517
    00518     GibbsSampler3(const vector<vector<double> >& pixel_intensities_,
    00519                  const vector<vector<vector<double> > >& spot_intensities_,
    00520                  const vector<Vector<4> >& spots_,
    00521                  const vector<vector<int> >& spot_pixels_,
    00522                  const Matrix<3> A_,
    00523                  const Vector<3> pi_,
    00524                  double variance_,
    00525                  int sample_iterations_)
    00526     :pixel_intensities(pixel_intensities_), //pixel_intensities: [frame][pixels]
    00527      spot_intensities(spot_intensities_),   //spot_intensities: [spot][frame][pixel]
    00528      spots(spots_),
    00529      spot_pixels(spot_pixels_),
    00530      A(A_),
    00531      pi(pi_),
    00532      base_variance(variance_),
    00533      variance(variance_),
    00534      sample_iterations(sample_iterations_),
    00535      num_frames(pixel_intensities.size()),
    00536      num_pixels(pixel_intensities[0].size()),
    00537         //Observations vector. As usual for this application, the observations are just
    00538         //numbered integers which refer to data held elsewhere.
    00539      O(sequence(num_frames)),
    00540         //Start all spots OFF, so the intensity is 0. OFF is 1 or 2, not 0!!!
    00541         //sample_list: [sample][spot][frame]: list of samples drawn using Gibbs sampling
    00542      current_sample(spots.size(), vector<State>(num_frames, 2)), //current sample [spot][frame]
    00543         //pixel intensities assosciated with the current sample [frame][pixel]
    00544      current_sample_intensities(num_frames, vector<double>(num_pixels))
    00545     {
    00546         //Check a bunch of stuff
    00547         assert_same_size(pixel_intensities);
    00548         assert_same_size(spot_intensities);
    00549     }
    00550
    00551     ///Update the noide variance. Used for adding thermal noise.
    00552     ///@param v noise variance.
    00553     void set_variance(double v)
    00554     {
    00555         variance = v;
    00556     }
    00557
    00558     ///Reset the gibbs sampler oro the initial state (all spots off)
    00559     void reset()
    00560     {
    00561         vector<State> off(num_frames, 2);
    00562         fill(current_sample.begin(), current_sample.end(), off);
    00563
    00564         vector<double> black(num_pixels);
    00565         fill(current_sample_intensities.begin(), current_sample_intensities.end(), black);
    00566         variance = base_variance;
    00567     }
    00568
    00569     ///Get the next sample
    00570     ///@param rng Random number generator
    00571     template<class T> void next(T& rng)
    00572     {
    00573
    00574         std::vector<array<double, 3> > delta3;
    00575         for(int j=0; j < sample_iterations; j++)
    00576             for(int k=0; k < (int) spots.size(); k++)
    00577             {
    00578                 //Subtract off the spot we're interested in.
    00579                 remove_spot(current_sample_intensities, spot_intensities[k], current_sample[k], spot_pixels[k]);
    00580
    00581                 //Now current_sample_intensities is the image value for every spot in every frame,
    00582                 //except the current spot, which is always set to off. This allows us to add it in
    00583                 //easily.
    00584
    00585                 SpotWithBackgroundDriftMasked B3(current_sample_intensities, spot_intensities[k], pixel_intensities, variance, spot_pixels[k]);
    00586                 forward_algorithm_delta2<3>(A, pi, B3, O, delta3);
    00587
    00588                 current_sample[k] = backward_sampling<3,State, T>(A, delta3, rng);
    00589                 //Put the newly sampled spot in
    00590                 add_spot(current_sample_intensities, spot_intensities[k], current_sample[k], spot_pixels[k]);
    00591             }
    00592     }
    00593 /*  void next()
    00594     {
    00595         RngDrand48 rng;
    00596         next(rng);
    00597     }
    00598 */
    00599     ///Retrieve the current sample
    00600     const vector<vector<State> >& sample() const
    00601     {
    00602         return current_sample;
    00603     }
    00604     ///Retrieve the intensities for the current sample
    00605     const vector<vector<double> >& sample_intensities() const
    00606     {
    00607         return current_sample_intensities;
    00608     }
    00609
    00610 };
    00611 #endif
    00612
    00613 }
    00614
    00615 using SampledMultispot::SpotWithBackground;
    00616 using SampledMultispot::remove_spot;
    00617 using SampledMultispot::add_spot;
    00618 using SampledMultispot::compute_spot_intensity;
    00619 using SampledMultispot::compute_spot_intensity_hessian;
    00620 using SampledMultispot::compute_spot_intensity_derivatives;
    00621 using SampledMultispot::sequence;
    00622 using SampledMultispot::GibbsSampler;
    00623 using SampledMultispot::GibbsSampler2;
    00624
    00625 #endif