00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
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
00043
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
00053
00054
00055
00056
00057
00058
00059
00060
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
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)
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)
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
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)
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)
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
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)
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)
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
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
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
00184
00185
00186
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
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
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_),
00234 spot_intensities(spot_intensities_),
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
00244
00245 O(sequence(num_frames)),
00246
00247
00248 current_sample(spots.size(), vector<State>(num_frames, 2)),
00249
00250 current_sample_intensities(num_frames, vector<double>(num_pixels))
00251 {
00252
00253 assert_same_size(pixel_intensities);
00254 assert_same_size(spot_intensities);
00255
00256 }
00257
00258
00259
00260 void set_variance(double v)
00261 {
00262 variance = v;
00263 }
00264
00265
00266
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
00278
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
00285 remove_spot(current_sample_intensities, spot_intensities[k], current_sample[k]);
00286
00287
00288
00289
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
00295 add_spot(current_sample_intensities, spot_intensities[k], current_sample[k]);
00296 }
00297 }
00298
00299
00300
00301
00302
00303
00304
00305 const vector<vector<State> >& sample() const
00306 {
00307 return current_sample;
00308 }
00309
00310 const vector<vector<double> >& sample_intensities() const
00311 {
00312 return current_sample_intensities;
00313 }
00314
00315 };
00316
00317
00318
00319
00320
00321
00322
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_),
00356 spot_intensities(spot_intensities_),
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
00367
00368 O(sequence(num_frames)),
00369
00370
00371 current_sample(spots.size(), vector<State>(num_frames, 2)),
00372
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
00378 assert_same_size(pixel_intensities);
00379 assert_same_size(spot_intensities);
00380 }
00381
00382
00383
00384 void set_variance(double v)
00385 {
00386 variance = v;
00387 }
00388
00389
00390
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
00402
00403 template<class T> void next(T& rng)
00404 {
00405
00406
00407
00408
00409
00410
00411
00412
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
00418
00419 remove_spot(current_sample_intensities, spot_intensities[k], current_sample[k], spot_pixels[k]);
00420
00421
00422
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436
00437
00438
00439
00440
00441
00442
00443
00444
00445
00446
00447
00448
00449
00450
00451
00452
00453 SpotWithBackgroundMasked B3(current_sample_intensities, spot_intensities[k], pixel_intensities, variance, spot_pixels[k]);
00454
00455 forward_algorithm_delta2<3>(A, pi, B3, O, delta3);
00456
00457
00458
00459
00460
00461
00462
00463
00464
00465
00466
00467
00468
00469 current_sample[k] = backward_sampling<3,State, T>(A, delta3, rng);
00470
00471
00472 add_spot(current_sample_intensities, spot_intensities[k], current_sample[k], spot_pixels[k]);
00473
00474 }
00475
00476 }
00477
00478
00479
00480
00481
00482
00483
00484 const vector<vector<State> >& sample() const
00485 {
00486 return current_sample;
00487 }
00488
00489 const vector<vector<double> >& sample_intensities() const
00490 {
00491 return current_sample_intensities;
00492 }
00493
00494 };
00495
00496 #if 0
00497
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_),
00527 spot_intensities(spot_intensities_),
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
00538
00539 O(sequence(num_frames)),
00540
00541
00542 current_sample(spots.size(), vector<State>(num_frames, 2)),
00543
00544 current_sample_intensities(num_frames, vector<double>(num_pixels))
00545 {
00546
00547 assert_same_size(pixel_intensities);
00548 assert_same_size(spot_intensities);
00549 }
00550
00551
00552
00553 void set_variance(double v)
00554 {
00555 variance = v;
00556 }
00557
00558
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
00570
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
00579 remove_spot(current_sample_intensities, spot_intensities[k], current_sample[k], spot_pixels[k]);
00580
00581
00582
00583
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
00590 add_spot(current_sample_intensities, spot_intensities[k], current_sample[k], spot_pixels[k]);
00591 }
00592 }
00593
00594
00595
00596
00597
00598
00599
00600 const vector<vector<State> >& sample() const
00601 {
00602 return current_sample;
00603 }
00604
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