ThreeB 1.1
|
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 #include <cstdlib> 00021 #include <cerrno> 00022 #include <cstring> 00023 #include <stack> 00024 #include <algorithm> 00025 #include <climits> 00026 #include <iomanip> 00027 #include <map> 00028 #include <tr1/memory> 00029 #include <cvd/image_io.h> 00030 #include <cvd/image_convert.h> 00031 #include <cvd/morphology.h> 00032 #include <cvd/connected_components.h> 00033 #include <cvd/draw.h> 00034 #include <cvd/vector_image_ref.h> 00035 00036 #include <cvd/random.h> 00037 #include <cvd/timer.h> 00038 #include <gvars3/instances.h> 00039 00040 #include <TooN/functions/derivatives.h> 00041 #include <TooN/determinant.h> 00042 #include <TooN/SymEigen.h> 00043 #include <TooN/optimization/conjugate_gradient.h> 00044 00045 #include "conjugate_gradient_only.h" 00046 #include "forward_algorithm.h" 00047 #include "numerical_derivatives.h" 00048 #include "storm.h" 00049 #include "storm_imagery.h" 00050 #include "debug.h" 00051 #include "sampled_multispot.h" 00052 #include "mt19937.h" 00053 #include "utility.h" 00054 #include "multispot5.h" 00055 00056 00057 //For benchmarking... 00058 #define TIME(X) 00059 //#define TIME(X) X 00060 00061 using namespace std; 00062 using namespace CVD; 00063 using namespace GVars3; 00064 using namespace TooN; 00065 using namespace std::tr1; 00066 00067 ///Empty destructor 00068 UserInterfaceCallback::~UserInterfaceCallback(){} 00069 00070 ///User interface callback class which does nothing. 00071 class NullUICallback: public UserInterfaceCallback 00072 { 00073 void per_spot(int, int, int, int){}; 00074 void per_modification(int, int, int){}; 00075 void per_pass(int, int, const std::vector<TooN::Vector<4> >&){}; 00076 void perhaps_stop(){}; 00077 }; 00078 00079 ///Factory function to generate an instance of NullGraphics 00080 ///@ingroup gStorm 00081 auto_ptr<UserInterfaceCallback> null_ui() 00082 { 00083 return auto_ptr<UserInterfaceCallback>(new NullUICallback); 00084 } 00085 00086 00087 //Declare the graphics classes. 00088 //These provide all the debug drawing operations so that the code can run in GUI and 00089 //headless mode easily and without macros. 00090 FitSpotsGraphics::~FitSpotsGraphics(){} 00091 00092 ///Graphics class which does absoloutely nothing 00093 ///@ingroup gStorm 00094 class NullGraphics: public FitSpotsGraphics 00095 { 00096 public: 00097 virtual void init(CVD::ImageRef){} 00098 virtual void draw_krap(const std::vector<TooN::Vector<4> >&, const CVD::Image<CVD::byte>&, const BBox&, int, TooN::Vector<4>){} 00099 virtual void swap(){} 00100 virtual void draw_pixels(const std::vector<CVD::ImageRef>&, float, float, float, float){} 00101 virtual void draw_bbox(const BBox&){} 00102 virtual void glDrawCross(const TooN::Vector<2>&, int){} 00103 virtual ~NullGraphics(){} 00104 }; 00105 00106 ///Factory function to generate an instance of NullGraphics 00107 ///@ingroup gStorm 00108 auto_ptr<FitSpotsGraphics> null_graphics() 00109 { 00110 return auto_ptr<FitSpotsGraphics>(new NullGraphics); 00111 } 00112 00113 00114 ///There are two sensible ways of storing the state vector of spot positions. 00115 ///This function converts between them. See also spots_to_vector. 00116 ///@param s list of spots to convert 00117 ///@ingroup gUtility 00118 Vector<> spots_to_Vector(const vector<Vector<4> >& s) 00119 { 00120 Vector<> r(s.size()*4); 00121 for(unsigned int i=0; i < s.size(); i++) 00122 { 00123 r[i*4+0] = s[i][0]; 00124 r[i*4+1] = s[i][1]; 00125 r[i*4+2] = s[i][2]; 00126 r[i*4+3] = s[i][3]; 00127 } 00128 return r; 00129 } 00130 00131 ///There are two sensible ways of storing the state vector of spot positions. 00132 ///This function converts between them. See also spots_to_Vector. 00133 ///@param s list of spots to convert 00134 ///@ingroup gUtility 00135 vector<Vector<4> > spots_to_vector(const Vector<>& s) 00136 { 00137 vector<Vector<4> > r(s.size()/4); 00138 for(unsigned int i=0; i < r.size(); i++) 00139 r[i] = s.slice<Dynamic, 4>(i*4, 4); 00140 return r; 00141 } 00142 00143 ///Normalize an image for display purposes. 00144 ///@ingroup gUtility 00145 Image<byte> scale_to_bytes(const Image<float>& im, float lo, float hi) 00146 { 00147 Image<byte> out(im.size()); 00148 for(int r=0; r < out.size().y-0; r++) 00149 for(int c=0; c < out.size().x-0; c++) 00150 out[r][c] = (int)floor((im[r][c]-lo)*255/(hi-lo)); 00151 return out; 00152 } 00153 00154 ///Normalize an image for display purposes. 00155 ///@ingroup gUtility 00156 Image<byte> scale_to_bytes(const Image<float>& im) 00157 { 00158 float lo = *min_element(im.begin(), im.end()); 00159 float hi = *max_element(im.begin(), im.end()); 00160 Image<byte> out(im.size()); 00161 for(int r=0; r < out.size().y-0; r++) 00162 for(int c=0; c < out.size().x-0; c++) 00163 out[r][c] = (int)floor((im[r][c]-lo)*255/(hi-lo)); 00164 00165 return out; 00166 } 00167 00168 ///Average the input image stack for display purposes 00169 ///@ingroup gUtility 00170 Image<float> average_image(const vector<Image<float> >& ims) 00171 { 00172 assert_same_size(ims); 00173 Image<float> r(ims[0].size(), 0); 00174 00175 for(unsigned int i=0; i < ims.size(); i++) 00176 transform(r.begin(), r.end(), ims[i].begin(), r.begin(), plus<float>()); 00177 00178 transform(r.begin(), r.end(), r.begin(), bind2nd(multiplies<float>(), 1./ims.size())); 00179 return r; 00180 } 00181 00182 00183 //////////////////////////////////////////////////////////////////////////////// 00184 //////////////////////////////////////////////////////////////////////////////// 00185 //////////////////////////////////////////////////////////////////////////////// 00186 //////////////////////////////////////////////////////////////////////////////// 00187 //////////////////////////////////////////////////////////////////////////////// 00188 //////////////////////////////////////////////////////////////////////////////// 00189 //////////////////////////////////////////////////////////////////////////////// 00190 //////////////////////////////////////////////////////////////////////////////// 00191 00192 ///Closure hoding the data required do use GibbsSampler2 00193 ///See FitSpots for naming of variables. 00194 ///@ingroup gStorm 00195 class DataForMCMC 00196 { 00197 protected: 00198 const vector<ImageRef>& pixels; 00199 const vector<vector<double> >& pixel_intensities;//[frame][pixel] 00200 const double mu_brightness, sigma_brightness, mu_blur, sigma_blur; 00201 const double variance; 00202 const int samples, sample_iterations; 00203 const Matrix<3> A; 00204 const Vector<3> pi; 00205 MT19937& rng; 00206 public: 00207 00208 MT19937& get_rng()const 00209 { 00210 return rng; 00211 } 00212 00213 public: 00214 DataForMCMC(const vector<ImageRef>& pixels_, 00215 const vector<vector<double> >& pixel_intensities_, 00216 double mu_brightness_, 00217 double sigma_brightness_, 00218 double mu_blur_, 00219 double sigma_blur_, 00220 double variance_, 00221 int samples_, 00222 int sample_iterations_, 00223 Matrix<3> A_, 00224 Vector<3> pi_, 00225 MT19937& rng_) 00226 :pixels(pixels_), 00227 pixel_intensities(pixel_intensities_), 00228 mu_brightness(mu_brightness_), 00229 sigma_brightness(sigma_brightness_), 00230 mu_blur(mu_blur_), 00231 sigma_blur(sigma_blur_), 00232 variance(variance_), 00233 samples(samples_), 00234 sample_iterations(sample_iterations_), 00235 A(A_), 00236 pi(pi_), 00237 rng(rng_) 00238 {} 00239 }; 00240 00241 ///Class implementing the Kahan summation algorithm 00242 ///to allow accurate summation of very large numbers of 00243 ///doubles. 00244 ///@ingroup gUtility 00245 class Kahan{ 00246 private: 00247 double y; ///< Input with the compensation removed. Temporary working space. 00248 double c; ///< Running compenstation for low-order bits. 00249 double t; ///< y + sum, which loses low order bits of y. Temporary working space. 00250 public: 00251 double sum; ///< running sum 00252 00253 Kahan() 00254 :c(0),sum(0) 00255 {} 00256 00257 /// Add a number to the running sum 00258 /// @param i Number to add 00259 void add(double i) 00260 { 00261 //y = input -c 00262 y = i; 00263 y-= c; 00264 00265 //t = sum + y 00266 t = sum; 00267 t += y; 00268 00269 //c = (t - sum) - y 00270 //c = ((sum + y) - sum) - y) 00271 c = t; 00272 c -= sum; 00273 c -= y; 00274 sum = t; 00275 } 00276 }; 00277 00278 00279 ///Class for computing the negitve free energy using thermodynamic integration. 00280 class NegativeFreeEnergy: public DataForMCMC 00281 { 00282 public: 00283 00284 ///@param d Necessary data 00285 NegativeFreeEnergy(const DataForMCMC& d) 00286 :DataForMCMC(d) 00287 { 00288 } 00289 00290 ///Give the noise variance given a sample number and growth parameters 00291 ///@param sample Sample number 00292 ///@param samples Total number of samples 00293 ///@param base_sigma Starting value of sigme 00294 ///@param scale_pow Exponent scaling 00295 double variance_from_sample(double sample, double samples, double base_sigma, double scale_pow) const 00296 { 00297 double scale = pow(1.25, sample * 1. / samples * scale_pow); 00298 double sigma = base_sigma * scale; 00299 double new_variance = sq(sigma); 00300 00301 return new_variance; 00302 } 00303 00304 ///Estimate free energy using the Slow Growth Thermodynamic Integration method given in 00305 ///Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford. M. Neal, 1993 00306 ///Except using a 5 point stencil instead of forward differenceing in Eq 6.30 00307 ///@param spots list of spots 00308 ///@param spot_pixels Mask around each spot to use to save on computation 00309 double compute_with_mask(const Vector<>& spots, const vector<vector<int> >& spot_pixels) const 00310 { 00311 //Estimate free energy using the Slow Growth Thermodynamic Integration method given in 00312 //Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford. M. Neal, 1993 00313 //Except using a 5 point stencil instead of forward differenceing in Eq 6.30 00314 double base_sigma = sqrt(variance); // should be 1 00315 double scale_pow = 100.0; 00316 00317 const unsigned int nspots = spots.size()/4; 00318 const unsigned int nframes = pixel_intensities.size(); 00319 const unsigned int npixels = pixels.size(); 00320 assert(spots.size() %4 == 0); 00321 assert(spot_pixels.size() == nspots); 00322 00323 //Compute the intensities and derivatives for all spot 00324 vector<vector<double> > spot_intensity; //[spot][pixel] 00325 for(unsigned int i=0; i < nspots; i++) 00326 spot_intensity.push_back(compute_spot_intensity(pixels, spots.slice<Dynamic,4>(i*4,4))); 00327 00328 GibbsSampler2 sampler(pixel_intensities, spot_intensity, spots_to_vector(spots), spot_pixels, A, pi, variance, sample_iterations); 00329 00330 double sum = 0; 00331 Kahan ksum; 00332 for(int sample=0; sample < samples; sample++) 00333 { 00334 //Compute the positions of the surrounding steps 00335 double var1 = variance_from_sample(sample-2, samples, base_sigma, scale_pow); 00336 double var2 = variance_from_sample(sample-1, samples, base_sigma, scale_pow); 00337 double var3 = variance_from_sample(sample+1, samples, base_sigma, scale_pow); 00338 double var4 = variance_from_sample(sample+2, samples, base_sigma, scale_pow); 00339 00340 //Take a sample 00341 sampler.set_variance(var2); 00342 sampler.next(DataForMCMC::get_rng()); 00343 00344 //Compute the SSD. This is fixed regardless of sigma. 00345 double err_sum=0; 00346 for(unsigned int frame=0; frame < nframes; frame++) 00347 for(unsigned int pixel=0; pixel < npixels; pixel++) 00348 err_sum -= sq(pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]); 00349 00350 //Compute the derivative using a five point stencil 00351 //This could be done in a better but less clear way 00352 double e1 = err_sum / (2*var1) - npixels*nframes*::log(2*M_PI*var1)/2; 00353 double e2 = err_sum / (2*var2) - npixels*nframes*::log(2*M_PI*var2)/2; 00354 double e3 = err_sum / (2*var3) - npixels*nframes*::log(2*M_PI*var3)/2; 00355 double e4 = err_sum / (2*var4) - npixels*nframes*::log(2*M_PI*var4)/2; 00356 sum += (-e1 + 8*e2 - 8*e3 + e4)/12; 00357 00358 ksum.add(-e1/12); 00359 ksum.add(8*e2/12); 00360 ksum.add(-8*e3/12); 00361 ksum.add(e4/12); 00362 } 00363 00364 double log_final = (log(variance_from_sample(samples, samples, base_sigma, scale_pow)*2*M_PI)/2) * npixels * nframes; 00365 00366 double priors=0; 00367 for(unsigned int i=0; i < nspots; i++) 00368 { 00369 priors += log_log_normal(spots[i*4+0], mu_brightness, sigma_brightness); 00370 priors += log_log_normal(spots[i*4+1], mu_blur, sigma_blur); 00371 } 00372 00373 /*cout << "Thermo:\n"; 00374 cout << "sum: " << sum -log_final << endl; 00375 cout << "kah: " << ksum.sum -log_final << endl; 00376 cout << "priors: " << priors + sum -log_final << endl; 00377 cout << " kah: " << priors + ksum.sum -log_final << endl; 00378 */ 00379 //cout << log_final << endl; 00380 //cout << sum + log_final << endl; 00381 00382 00383 sampler.set_variance(variance); 00384 return -(sum+priors - log_final); 00385 } 00386 00387 ///Estimate free energy using the Slow Growth Thermodynamic Integration method given in 00388 ///Probalistic Inference Using Markov Chain Monte Carlo Methods, Radford. M. Neal, 1993 00389 ///Except using a 5 point stencil instead of forward differenceing in Eq 6.30 00390 ///@param spots list of spots 00391 double operator()(const Vector<>& spots) const 00392 { 00393 double base_sigma = sqrt(variance); // should be 1 00394 double scale_pow = 100.0; 00395 00396 const unsigned int nspots = spots.size()/4; 00397 const unsigned int nframes = pixel_intensities.size(); 00398 const unsigned int npixels = pixels.size(); 00399 assert(spots.size() %4 == 0); 00400 00401 //Compute the intensities and derivatives for all spot 00402 vector<vector<double> > spot_intensity; //[spot][pixel] 00403 for(unsigned int i=0; i < nspots; i++) 00404 spot_intensity.push_back(compute_spot_intensity(pixels, spots.slice<Dynamic,4>(i*4,4))); 00405 00406 GibbsSampler sampler(pixel_intensities, spot_intensity, spots_to_vector(spots), A, pi, variance, sample_iterations); 00407 00408 double sum = 0; 00409 Kahan ksum; 00410 for(int sample=0; sample < samples; sample++) 00411 { 00412 //Compute the positions of the surrounding steps 00413 double var1 = variance_from_sample(sample-2, samples, base_sigma, scale_pow); 00414 double var2 = variance_from_sample(sample-1, samples, base_sigma, scale_pow); 00415 double var3 = variance_from_sample(sample+1, samples, base_sigma, scale_pow); 00416 double var4 = variance_from_sample(sample+2, samples, base_sigma, scale_pow); 00417 00418 //Take a sample 00419 sampler.set_variance(var2); 00420 sampler.next(DataForMCMC::get_rng()); 00421 00422 //Compute the SSD. This is fixed regardless of sigma. 00423 double err_sum=0; 00424 for(unsigned int frame=0; frame < nframes; frame++) 00425 for(unsigned int pixel=0; pixel < npixels; pixel++) 00426 err_sum -= sq(pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]); 00427 00428 //Compute the derivative using a five point stencil 00429 //This could be done in a better but less clear way 00430 double e1 = err_sum / (2*var1) - npixels*nframes*::log(2*M_PI*var1)/2; 00431 double e2 = err_sum / (2*var2) - npixels*nframes*::log(2*M_PI*var2)/2; 00432 double e3 = err_sum / (2*var3) - npixels*nframes*::log(2*M_PI*var3)/2; 00433 double e4 = err_sum / (2*var4) - npixels*nframes*::log(2*M_PI*var4)/2; 00434 sum += (-e1 + 8*e2 - 8*e3 + e4)/12; 00435 00436 ksum.add(-e1/12); 00437 ksum.add(8*e2/12); 00438 ksum.add(-8*e3/12); 00439 ksum.add(e4/12); 00440 } 00441 00442 double log_final = (log(variance_from_sample(samples, samples, base_sigma, scale_pow)*2*M_PI)/2) * npixels * nframes; 00443 00444 double priors=0; 00445 for(unsigned int i=0; i < nspots; i++) 00446 { 00447 priors += log_log_normal(spots[i*4+0], mu_brightness, sigma_brightness); 00448 priors += log_log_normal(spots[i*4+1], mu_blur, sigma_blur); 00449 } 00450 00451 /*cout << "Thermo:\n"; 00452 cout << "sum: " << sum -log_final << endl; 00453 cout << "kah: " << ksum.sum -log_final << endl; 00454 cout << "priors: " << priors + sum -log_final << endl; 00455 cout << " kah: " << priors + ksum.sum -log_final << endl; 00456 */ 00457 //cout << log_final << endl; 00458 //cout << sum + log_final << endl; 00459 00460 00461 sampler.set_variance(variance); 00462 return -(sum+priors - log_final); 00463 } 00464 }; 00465 00466 /// Class for sorting a list of indexes to an array of spots lexicographically 00467 /// according to the 2D positions of the spots. 00468 /// 00469 ///@param Cmp comparator function to specify less or greater 00470 ///@param First most significant position index 00471 ///@ingroup gMultiSpot 00472 template<class Cmp, int First> struct IndexLexicographicPosition{ 00473 const vector<Vector<4> >& spots; ///< Keep around the array of spots for later comprison 00474 00475 /// @param s Vector to sort indices of 00476 IndexLexicographicPosition(const vector<Vector<4> >& s) 00477 :spots(s) 00478 {} 00479 00480 00481 static const int Second = First==2?3:2; ///< Second most siginifcant position index for sorting 00482 00483 /// Compare two indexes into the array of spots 00484 bool operator()(int a, int b) 00485 { 00486 Cmp cmp; 00487 00488 if(cmp(spots[a][First], spots[b][First])) 00489 return true; 00490 else if(spots[a][First] == spots[b][First]) 00491 return cmp(spots[a][Second], spots[b][Second]); 00492 else 00493 return false; 00494 } 00495 }; 00496 00497 00498 ///Closure holding image data generated using samples drawn from the model. 00499 ///NB this is used with one spot removed (i.e. set to dark). As a result, the data 00500 ///is treated as the background when considering that spot in isolation. 00501 ///See FitSpots for naming of variables. 00502 ///@ingroup gStorm 00503 struct SampledBackgroundData 00504 { 00505 const vector<vector<vector<double> > >& sample_intensities_without_spot; //[sample][frame][pixel] 00506 const vector<vector<double> >& pixel_intensities; //[frame][pixel] 00507 const vector<ImageRef> pixels; 00508 00509 double mu_brightness, sigma_brightness, mu_blur, sigma_blur; 00510 const Matrix<3> A; 00511 const Vector<3> pi; 00512 double variance; 00513 00514 const vector<int> O; 00515 00516 SampledBackgroundData( 00517 const vector<vector<vector<double> > >& sample_intensities_without_spot_, 00518 const vector<vector<double> >& pixel_intensities_, 00519 const vector<ImageRef> pixels_, 00520 double mu_brightness_, 00521 double sigma_brightness_, 00522 double mu_blur_, 00523 double sigma_blur_, 00524 const Matrix<3> A_, 00525 const Vector<3> pi_, 00526 double variance_) 00527 :sample_intensities_without_spot(sample_intensities_without_spot_), 00528 pixel_intensities(pixel_intensities_), 00529 pixels(pixels_), 00530 mu_brightness(mu_brightness_), 00531 sigma_brightness(sigma_brightness_), 00532 mu_blur(mu_blur_), 00533 sigma_blur(sigma_blur_), 00534 A(A_), 00535 pi(pi_), 00536 variance(variance_), 00537 O(sequence(pixel_intensities.size())) 00538 { 00539 } 00540 }; 00541 00542 //!@cond Doxygen_Suppress 00543 //Do not use. 00544 struct SpotProbabilityWithSampledBackgroundFAKE: public SampledBackgroundData 00545 { 00546 SpotProbabilityWithSampledBackgroundFAKE(const SampledBackgroundData& d) 00547 :SampledBackgroundData(d) 00548 { 00549 } 00550 00551 //Compute the probability of spot 00552 double operator()(const Vector<4>& spot) const 00553 { 00554 vector<double> spot_intensities = compute_spot_intensity(pixels, spot); 00555 00556 double sum_log_prob = 0; 00557 00558 for(unsigned int s=0; s < sample_intensities_without_spot.size(); s++) 00559 { 00560 SpotWithBackground B(sample_intensities_without_spot[s], spot_intensities, pixel_intensities, variance); 00561 double log_prob = forward_algorithm(A, pi, B, O); 00562 00563 double logprior = log_log_normal(spot[0], mu_brightness, sigma_brightness) + 00564 log_log_normal(spot[1], mu_blur, sigma_blur); 00565 00566 //cout << setprecision(10) <<fixed << "Prob : " << log_prob + logprior << endl; 00567 00568 sum_log_prob += log_prob + logprior; 00569 } 00570 00571 double average_log_prob = sum_log_prob / sample_intensities_without_spot.size(); 00572 00573 // cout << "Ave Prob : " << average_log_prob << endl; 00574 00575 if(spot[0] < 0 || spot[1] < 0) 00576 return 1e100; 00577 00578 return average_log_prob; 00579 } 00580 }; 00581 //@endcond 00582 00583 ///Compute the derivative of the negative log probability with respect 00584 ///to the parameters of one spot, given some samples of the other spots. 00585 ///@ingroup gStorm 00586 struct SpotNegProbabilityDiffWithSampledBackground: public SampledBackgroundData 00587 { 00588 ///@param d Necessary data for construction 00589 SpotNegProbabilityDiffWithSampledBackground(const SampledBackgroundData& d) 00590 :SampledBackgroundData(d) 00591 { 00592 } 00593 00594 ///Compute the probability of spot 00595 ///@param spot Spot position 00596 Vector<4> operator()(const Vector<4>& spot) const 00597 { 00598 if(spot[0] <= 0 || spot[1] <= 0) 00599 return Ones * std::numeric_limits<double>::quiet_NaN(); 00600 00601 vector<pair<double, Vector<4> > > spot_intensities = compute_spot_intensity_derivatives(pixels, spot); 00602 00603 Vector<4> sum_diff_log = Zeros; 00604 00605 for(unsigned int s=0; s < sample_intensities_without_spot.size(); s++) 00606 { 00607 SpotWithBackground B(sample_intensities_without_spot[s], spot_intensities, pixel_intensities, variance); 00608 00609 pair<double, Vector<4> > r = forward_algorithm_deriv(A, pi, B, O); 00610 00611 sum_diff_log += r.second; 00612 } 00613 00614 Vector<4> diff_log = sum_diff_log / sample_intensities_without_spot.size(); 00615 00616 //Compute the log probability of the prior 00617 Vector<4> logprior_deriv = makeVector(diff_log_log_normal(spot[0], mu_brightness, sigma_brightness), 00618 diff_log_log_normal(spot[1], mu_blur, sigma_blur), 0, 0); 00619 00620 return -(diff_log + logprior_deriv); 00621 } 00622 }; 00623 00624 00625 ///Class for computing the Hessian of the negative free energy. 00626 ///@ingroup gStorm 00627 class FreeEnergyHessian: public DataForMCMC 00628 { 00629 public: 00630 00631 ///Constructor 00632 ///@param d All data required 00633 FreeEnergyHessian(const DataForMCMC& d) 00634 :DataForMCMC(d) 00635 { 00636 } 00637 00638 ///Compute the Hessian 00639 ///@param spots All spot positions 00640 ///@param spot spot to compute hessian for 00641 Matrix<4> hessian(const vector<Vector<4> >& spots, int spot) const 00642 { 00643 cout << "----Computing pure MCMC hessian\n"; 00644 const unsigned int nspots = spots.size(); 00645 const unsigned int nframes = pixel_intensities.size(); 00646 const unsigned int npixels = pixels.size(); 00647 cout << spot << " " << nspots << " " << nframes << " " << npixels << endl; 00648 00649 vector<vector<double> > spot_intensity; //[spot][pixel] 00650 for(unsigned int i=0; i < nspots; i++) 00651 spot_intensity.push_back(compute_spot_intensity(pixels, spots[i])); 00652 00653 vector<tuple<double, Vector<4>, Matrix<4> > > spot_hess_etc = compute_spot_intensity_hessian(pixels, spots[spot]); 00654 00655 GibbsSampler sampler(pixel_intensities, spot_intensity, spots, A, pi, variance, sample_iterations); 00656 00657 //Compute derivative of log probability, summed (ie integrated) 00658 Matrix<4> sum_hess1 = Zeros(spots.size()); 00659 Matrix<4> sum_hess2 = Zeros(spots.size()); 00660 Vector<4> sum_deriv = Zeros(spots.size()); 00661 00662 for(int sample=0; sample < samples; sample++) 00663 { 00664 sampler.next(rng); 00665 00666 //Compute d log P(data | x, model) / d model, for a given sample 00667 //And the hessian 00668 Matrix<4> hess = Zeros(spots.size()); 00669 Vector<4> deriv = Zeros(spots.size()); 00670 for(unsigned int frame=0; frame < nframes; frame++) 00671 { 00672 for(unsigned int pixel=0; pixel < npixels; pixel++) 00673 { 00674 double e = pixel_intensities[frame][pixel] - sampler.sample_intensities()[frame][pixel]; 00675 //Build up the derivative 00676 if(sampler.sample()[spot][frame] == 0) 00677 { 00678 hess += e * get<2>(spot_hess_etc[pixel]) - get<1>(spot_hess_etc[pixel]).as_col() * get<1>(spot_hess_etc[pixel]).as_row(); 00679 deriv += e * get<1>(spot_hess_etc[pixel]); 00680 00681 } 00682 } 00683 } 00684 00685 hess[0][0] += hess_log_log_normal(spots[spot][0], mu_brightness, sigma_brightness); 00686 hess[1][1] += hess_log_log_normal(spots[spot][1], mu_blur, sigma_blur); 00687 sum_hess1 += hess; 00688 00689 deriv[0] += diff_log_log_normal(spots[spot][0], mu_brightness, sigma_brightness); 00690 deriv[1] += diff_log_log_normal(spots[spot][1], mu_blur, sigma_blur); 00691 00692 sum_hess2 += deriv.as_col() * deriv.as_row(); 00693 00694 sum_deriv = sum_deriv + deriv; 00695 } 00696 00697 sum_hess1 /= (samples * variance); 00698 sum_hess2 /= (samples * variance); 00699 sum_deriv /= (samples * variance); 00700 00701 00702 cout << sum_hess1 << endl; 00703 cout << sum_hess2 << endl; 00704 cout << sum_deriv.as_col() * sum_deriv.as_row() << endl; 00705 00706 cout << "......." << sum_deriv << endl; 00707 //Put in the prior 00708 //The derivative prior parts cancel out. 00709 //Rather sensibly this means that the second derivatives can be 00710 //computed without reference to the prior, and then the prior can 00711 //be added in later, i.e.: P(data, parameters) = P(data|parameter) P(parameters) 00712 00713 //The second derivatives have been constructed to be diagonal 00714 DiagonalMatrix<4> hess_prior = Zeros(spots.size()); 00715 00716 cout << "sum of parts = \n" << sum_hess1 + sum_hess2 - sum_deriv.as_col() * sum_deriv.as_row() << endl; 00717 //TooN cannot currently add DiagonalMatrix to Matrix!! 00718 //sum_hess.diagonal_slice() += hess_prior.diagonal_slice(); 00719 00720 cout << "++++Done Computing pure MCMC hessian\n"; 00721 return sum_hess1 + sum_hess2 - sum_deriv.as_col() * sum_deriv.as_row(); 00722 } 00723 }; 00724 00725 ///Comparator functor for the first element of a std::pair 00726 struct LessSecond 00727 { 00728 ///Comparison function 00729 ///@param a 00730 ///@param b 00731 template<class A, class B> bool operator()(const pair<A,B>& a, const pair<A,B>& b) const 00732 { 00733 return a.second < b.second; 00734 } 00735 }; 00736 00737 //!@cond Doxygen_Suppress 00738 struct NthDeriv{ 00739 00740 const SpotNegProbabilityDiffWithSampledBackground& compute_deriv; 00741 int i; 00742 00743 NthDeriv(const SpotNegProbabilityDiffWithSampledBackground& c, int ii) 00744 :compute_deriv(c),i(ii) 00745 {} 00746 00747 double operator()(const Vector<4>& f) const 00748 { 00749 return compute_deriv(f)[i]; 00750 } 00751 }; 00752 //!@endcond 00753 00754 ///Compute the Hessian of the log probability. The background is sampled rather sparsely, and the 00755 ///spot in question is sampled much more densely using FFBS. 00756 ///@param spot Spot parameters 00757 ///@param d Background brightness (from other spots) 00758 ///@param bs_iterations Exter backward sampling iterations 00759 ///@param rng Random number generator 00760 ///@returns the Hessian of the log probability around the spot 00761 Matrix<4> sampled_background_spot_hessian_ffbs(const Vector<4>& spot, const SampledBackgroundData& d, int bs_iterations, MT19937& rng) 00762 { 00763 vector<tuple<double, Vector<4>, Matrix<4> > > spot_hess_etc = compute_spot_intensity_hessian(d.pixels, spot); 00764 vector<double> spot_intensities = compute_spot_intensity(d.pixels, spot); 00765 00766 Matrix<4> sum_hess_log = Zeros; 00767 Matrix<4> sum_diff2_log = Zeros; 00768 00769 vector<State> current_sample; 00770 00771 const unsigned int nframes = d.pixel_intensities.size(); 00772 const unsigned int npixels = d.pixels.size(); 00773 00774 Matrix<4> sum_hess = Zeros; 00775 Vector<4> sum_deriv = Zeros; 00776 00777 vector<pair<Matrix<4>, Vector<4> > > hess_and_deriv_part(nframes); 00778 00779 for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++) 00780 { 00781 SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensities, d.pixel_intensities, d.variance); 00782 00783 //Compute what the per-frame hess and deriv parts are 00784 //if the spot is on in a frame. 00785 for(unsigned int frame=0; frame < nframes; frame++) 00786 { 00787 Matrix<4> hess = Zeros; 00788 Vector<4> deriv = Zeros; 00789 00790 for(unsigned int pixel=0; pixel < npixels; pixel++) 00791 { 00792 double e = d.pixel_intensities[frame][pixel] - (d.sample_intensities_without_spot[s][frame][pixel] + spot_intensities[pixel]); 00793 //Build up the derivative 00794 hess += e * get<2>(spot_hess_etc[pixel]) - get<1>(spot_hess_etc[pixel]).as_col() * get<1>(spot_hess_etc[pixel]).as_row(); 00795 deriv += e * get<1>(spot_hess_etc[pixel]); 00796 } 00797 hess_and_deriv_part[frame] = make_pair(hess, deriv); 00798 } 00799 00800 //Forward filtering 00801 std::vector<array<double, 3> > delta = forward_algorithm_delta(d.A, d.pi, B, d.O); 00802 00803 for(int i=0; i < bs_iterations; i++) 00804 { 00805 current_sample = backward_sampling<3,State>(d.A, delta, rng); 00806 00807 Matrix<4> hess = Zeros; 00808 Vector<4> deriv = Zeros; 00809 for(unsigned int frame=0; frame < nframes; frame++) 00810 if(current_sample[frame] == 0) 00811 { 00812 hess += hess_and_deriv_part[frame].first; 00813 deriv += hess_and_deriv_part[frame].second; 00814 } 00815 00816 sum_hess += hess + deriv.as_col() * deriv.as_row(); 00817 sum_deriv += deriv; 00818 } 00819 } 00820 00821 sum_hess /= (bs_iterations * d.sample_intensities_without_spot.size() * d.variance); 00822 sum_deriv /= (bs_iterations * d.sample_intensities_without_spot.size() * d.variance); 00823 00824 sum_hess -= sum_deriv.as_col() * sum_deriv.as_row(); 00825 00826 sum_hess[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness); 00827 sum_hess[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur); 00828 sum_deriv[0] += diff_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness); 00829 sum_deriv[1] += diff_log_log_normal(spot[1], d.mu_blur, d.sigma_blur); 00830 00831 //cout << "Turboderiv:" << sum_deriv << endl; 00832 //cout << "Turbohess:\n" << sum_hess << endl; 00833 00834 return sum_hess; 00835 } 00836 00837 ///Debugging function. Not mathematically correct. Do not use. 00838 Matrix<4> sampled_background_spot_hessian2(const Vector<4>& spot, const SampledBackgroundData& d) 00839 { 00840 vector<tuple<double, Vector<4>, Matrix<4> > > spot_intensities = compute_spot_intensity_hessian(d.pixels, spot); 00841 00842 Matrix<4> sum_hess_log = Zeros; 00843 Matrix<4> sum_diff2_log = Zeros; 00844 00845 for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++) 00846 { 00847 SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensities, d.pixel_intensities, d.variance); 00848 00849 double prob; 00850 Vector<4> diff; 00851 Matrix<4> hess; 00852 00853 tie(prob, diff, hess) = forward_algorithm_hessian(d.A, d.pi, B, d.O); 00854 00855 sum_hess_log += hess; 00856 00857 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); 00858 sum_diff2_log += diff.as_col() * diff.as_row(); 00859 } 00860 00861 Matrix<4> hess_log = sum_hess_log / d.sample_intensities_without_spot.size(); 00862 Matrix<4> diff2_log = sum_diff2_log / d.sample_intensities_without_spot.size(); 00863 00864 //Add in the prior 00865 00866 hess_log[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness); 00867 hess_log[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur); 00868 00869 return hess_log + diff2_log; 00870 } 00871 00872 ///Debugging function. Not mathematically correct. Do not use. 00873 Matrix<4> sampled_background_spot_hessian_FAKE(const Vector<4>& spot, const SampledBackgroundData& d) 00874 { 00875 vector<tuple<double, Vector<4>, Matrix<4> > > spot_intensities = compute_spot_intensity_hessian(d.pixels, spot); 00876 00877 Matrix<4> sum_hess_log = Zeros; 00878 00879 for(unsigned int s=0; s < d.sample_intensities_without_spot.size(); s++) 00880 { 00881 SpotWithBackground B(d.sample_intensities_without_spot[s], spot_intensities, d.pixel_intensities, d.variance); 00882 00883 double prob; 00884 Vector<4> diff; 00885 Matrix<4> hess; 00886 00887 tie(prob, diff, hess) = forward_algorithm_hessian(d.A, d.pi, B, d.O); 00888 00889 sum_hess_log += hess; 00890 } 00891 00892 Matrix<4> hess_log = sum_hess_log / d.sample_intensities_without_spot.size(); 00893 00894 //Add in the prior 00895 00896 hess_log[0][0] += hess_log_log_normal(spot[0], d.mu_brightness, d.sigma_brightness); 00897 hess_log[1][1] += hess_log_log_normal(spot[1], d.mu_blur, d.sigma_blur); 00898 00899 return hess_log; 00900 } 00901 00902 ///Which pixels belong to a given spot? 00903 ///Find the indices of those pixels 00904 ///@ingroup gStorm 00905 void get_spot_pixels(const vector<ImageRef>& pixels, const Vector<4>& spot, vector<int>& out) 00906 { 00907 //Go out to three sigma 00908 00909 vector<ImageRef> pix = getDisc(spot[1]*6 + 1); 00910 out.resize(0); 00911 ImageRef offset = ir_rounded(spot.slice<2,2>()); 00912 for(unsigned int j=0; j < pix.size(); j++) 00913 { 00914 int pos = lower_bound(pixels.begin(), pixels.end(), pix[j] + offset) - pixels.begin(); 00915 if(pos != (int)pixels.size() && pixels[pos] == pix[j] + offset) 00916 out.push_back(pos); 00917 } 00918 00919 if(out.size() == 0) 00920 { 00921 cout << "********************************\n"; 00922 cout << "********************************\n"; 00923 cout << "********************************\n"; 00924 cout << "********************************\n"; 00925 cout << "********************************\n"; 00926 cout << "Oe noes!11one\n"; 00927 cout << pix.size() << endl; 00928 } 00929 } 00930 00931 00932 ///Tokenize a line 00933 ///@ingroup gUtility 00934 vector<string> split(const string& line) 00935 { 00936 vector<string> v; 00937 istringstream i(line); 00938 string s; 00939 00940 while(!i.eof()) 00941 { 00942 i >> s; 00943 if(i.fail()) 00944 break; 00945 v.push_back(s); 00946 } 00947 return v; 00948 } 00949 00950 ///Generic version of itoa. 00951 ///How many times has this been reimplemented?? 00952 ///@param x Value to convert to string 00953 ///@ingroup gUtility 00954 template<class C> inline string xtoa(const C& x) 00955 { 00956 ostringstream os; 00957 os << x; 00958 return os.str(); 00959 } 00960 00961 ///Inverse of xtoa() 00962 ///How many times has this been reimplemented?? 00963 ///@param s String to convert 00964 ///@param msg Mesage to print on error 00965 ///@ingroup gUtility 00966 template<class C> inline C atox(const string& s, const string& msg) 00967 { 00968 C c; 00969 istringstream i(s); 00970 i >> c; 00971 if(i.fail()) 00972 throw LogFileParseError("Error parsing " + msg + ". Text is `" + s + "'."); 00973 return c; 00974 } 00975 00976 /**Parser for multispot 5 log files. 00977 00978 Log files are mostly line oriented and contain various records 00979 00980 The main records are: 00981 00982 Iteraton: \#ITNUM 00983 MAIN: <list of spot parameters> 00984 00985 Pass: \#PASSNUM 00986 MT19937 <random number generator state> 00987 PASS\#PASSNUM: <list of spot parameters> 00988 ENDCHECKPOINT 00989 00990 Note that MAIN is redundant since it will be the same as the following PASS 1 (or the first pass 00991 computed if restoring from a checkpoint). 00992 00993 Data should only be considered valid after ENDCHECKPOINT has been read 00994 00995 Iteration is written once per iteration, not once per pass. (FIXME) 00996 00997 Which moron invented this file format??? 00998 00999 Note that the file format hasn't beren fixed, so that the output can easily be compared to the 01000 output of the historic version which is known to be good. 01001 01002 @param in Stream to parse file from 01003 @ingroup gStorm 01004 */ 01005 StateParameters parse_log_file(istream& in) 01006 { 01007 //A line read from the file 01008 string line; 01009 01010 //State lines known to be OK 01011 string rngline, passline, iterationline; 01012 bool state_ok=0; 01013 01014 //State lines read in, with flags of goodness 01015 string new_rngline, new_passline, new_iterationline; 01016 bool new_rngline_ok=0, new_passline_ok=0, new_iterationline_ok=0; 01017 01018 unsigned int lineno=0; 01019 bool doing_gvars = 0; 01020 01021 vector<ImageRef> pixels; 01022 01023 while(!in.eof()) 01024 { 01025 getline(in, line); 01026 if(in.fail()) 01027 break; 01028 01029 lineno++; 01030 01031 if(line == "ENDGVARLIST") 01032 { 01033 if(!doing_gvars) 01034 throw LogFileParseError("Spurious end of GVars"); 01035 doing_gvars = 0; 01036 } 01037 else if(doing_gvars) 01038 { 01039 GUI.ParseLine(line); 01040 } 01041 else if(line == "BEGINGVARLIST") 01042 { 01043 doing_gvars = 1; 01044 } 01045 if(line.substr(0, 11) == "Iteration: ") 01046 { 01047 new_iterationline = line; 01048 new_iterationline_ok = true; 01049 } 01050 else if(line.substr(0, 4) == "PASS") 01051 { 01052 new_passline = line; 01053 if(new_passline_ok) 01054 throw LogFileParseError("Duplicate PASS on line " + xtoa(lineno)); 01055 new_passline_ok = true; 01056 } 01057 else if(line.substr(0, 8) == "MT19937 ") 01058 { 01059 new_rngline = line; 01060 if(new_rngline_ok) 01061 throw LogFileParseError("Duplicate MT19937 on line " + xtoa(lineno)); 01062 01063 new_rngline_ok = true; 01064 01065 } 01066 else if(line == "ENDCHECKPOINT") 01067 { 01068 if(new_passline_ok && new_rngline_ok && new_iterationline_ok) 01069 { 01070 iterationline = new_iterationline; 01071 rngline = new_rngline; 01072 passline = new_passline; 01073 } 01074 else 01075 throw LogFileParseError("Reached checkpoint with missing data: " 01076 "it=" + xtoa(new_iterationline_ok) + 01077 " pa=" + xtoa(new_passline_ok) + 01078 " rg=" + xtoa(new_rngline_ok) + " on line " + xtoa(lineno)); 01079 01080 //Don't reset iteration since it only appears once for each entire 01081 //set of passes. 01082 new_rngline_ok = 0; 01083 new_passline_ok = 0; 01084 01085 state_ok = true; 01086 } 01087 else if(line.substr(0, 7) == "PIXELS ") 01088 { 01089 vector<string> l = split(line); 01090 if( (l.size() - 1)%2 == 0) 01091 { 01092 int n = (l.size()-1)/2; 01093 pixels.resize(n); 01094 for(int i=0; i < n; i++) 01095 { 01096 pixels[i].x = atox<int>(l[i*2+1+0], "pixels"); 01097 pixels[i].y = atox<int>(l[i*2+1+1], "pixels"); 01098 } 01099 } 01100 else 01101 throw LogFileParseError("Bad PIXELS line"); 01102 } 01103 } 01104 01105 if(!state_ok) 01106 throw LogFileParseError("No state found"); 01107 01108 if(pixels.size() == 0) 01109 throw LogFileParseError("No pixels, or pixels is empty"); 01110 01111 //Now parse the lines 01112 StateParameters p; 01113 vector<string> l; 01114 01115 //Parse the iterations 01116 l = split(iterationline); 01117 p.iteration = atox<int>(l[1], "iteration"); 01118 01119 //Parse the random number generator 01120 p.rng =shared_ptr<MT19937>(new MT19937); 01121 { 01122 istringstream rng_s(rngline); 01123 try{ 01124 p.rng->read(rng_s); 01125 } 01126 catch(MT19937::ParseError p) 01127 { 01128 throw LogFileParseError("Error parsing MT19937"); 01129 } 01130 } 01131 01132 //Parse PASS and the listing of spots 01133 l = split(passline); 01134 if( (l.size() - 1 ) % 4 == 0) 01135 { 01136 p.pass = atox<int>(l[0].substr(4), "pass"); 01137 01138 for(unsigned int i=0; i < (l.size()-1)/4; i++) 01139 { 01140 cout << l[i*4+1+0] << endl; 01141 cout << l[i*4+1+1] << endl; 01142 cout << l[i*4+1+2] << endl; 01143 cout << l[i*4+1+3] << endl; 01144 p.spots.push_back(makeVector( 01145 atox<double>(l[i*4+1+0], "spot"), 01146 atox<double>(l[i*4+1+1], "spot"), 01147 atox<double>(l[i*4+1+2], "spot"), 01148 atox<double>(l[i*4+1+3], "spot"))); 01149 } 01150 01151 } 01152 else 01153 throw LogFileParseError("Wrong number of elements in PASS line"); 01154 01155 //Set up the pixels (oh god the pixels) 01156 p.pixels = pixels; 01157 01158 return p; 01159 } 01160 01161 01162 01163 ///Setup the parameters for a run using the old and deeply peculiar method. 01164 ///This includes the unpleasant and diffucult so use de-checkpointing code. 01165 ///wtf. The use of this function is very strongly deprecated. 01166 ///@param log_ratios Image from which region is selected. 01167 ///@param ims Input data 01168 ///@param pixels Region for spot fitting to run in 01169 StateParameters generate_state_parameters_ye_olde(const BasicImage<double>& log_ratios, const vector<Image<float> >& ims, vector<ImageRef> pixels){ 01170 sort(pixels.begin(), pixels.end()); 01171 01172 const double variance = 1; // it should be 01173 01174 //To scale the X axis of a log-normal distribution, only 01175 //the mu parameter needs to be changed... 01176 const double intensity_mu = GV3::get<double>("intensity.rel_mu", 0., -1) + log(sqrt(variance)); 01177 const double intensity_sigma = GV3::get<double>("intensity.rel_sigma", 0., -1); 01178 const double blur_mu = GV3::get<double>("blur.mu", 0., -1); 01179 const double blur_sigma = GV3::get<double>("blur.sigma", 0., -1); 01180 01181 //The region was extracted at a certain threshold. 01182 //These regions may be too small, so some post-region finding dilation 01183 //may be performed. New spots are only placed at pixels which exceed the threshold. 01184 //post_dilate.threshold is (if set) used as the placing threshold so the placing threshold 01185 //can be different from the region-finding threshold. 01186 01187 //Note that as a result of dliation, regions of <pixels> may be below the threshold. 01188 //In the historic version, this could affect new spot placement. This feature is not supported 01189 //in this version. 01190 double threshold = GV3::get<double>("threshold", 0, -1); 01191 const double post_threshold = GV3::get<double>("post_dilate.threshold", -1, 1); 01192 if(post_threshold != -1) 01193 threshold = post_threshold; 01194 01195 01196 //If dilation after region finding is to be performed, then do it here. 01197 const double post_dilate_radius = GV3::get<double>("post_dilate.radius", 0, -1); 01198 if(post_dilate_radius != 0) 01199 { 01200 Image<byte> pix(ims[0].size()); 01201 pix.fill(0); 01202 01203 for(unsigned int i=0; i < pixels.size(); i++) 01204 pix[pixels[i]] = 255; 01205 01206 Image<byte> dilated = morphology(pix, getDisc(post_dilate_radius), Morphology::BinaryDilate<byte>()); 01207 01208 pixels.clear(); 01209 01210 ImageRef p(0,0); 01211 do 01212 if(dilated[p]) 01213 pixels.push_back(p); 01214 while(p.next(dilated.size())); 01215 } 01216 01217 01218 assert_same_size(ims); 01219 if(log_ratios.size() != ims[0].size()) 01220 { 01221 cerr << "Bad log ratios size\n"; 01222 exit(1); 01223 } 01224 01225 vector<Vector<4> > spots; 01226 //Spots can be either put down automatically, or specified 01227 //The auto-initialization is very strange. 01228 if(GV3::get<bool>("spots.auto_initialise", 1, 1)) 01229 { 01230 //You never get two spots in the same disc in the second stage of the algorithm 01231 vector<ImageRef> disc = getDisc(GV3::get<double>("spot_spread", 3.1, 1)); 01232 01233 //Record all the pixels 01234 map<ImageRef, double> valid_pixels; 01235 for(unsigned int i=0; i < pixels.size(); i++) 01236 if(log_ratios[pixels[i]] > threshold) 01237 valid_pixels.insert(make_pair(pixels[i], log_ratios[pixels[i]])); 01238 01239 01240 //Get some initial spots by finding the local maxima 01241 ImageRef neighbours[8] = { 01242 ImageRef(-1, -1), 01243 ImageRef( 0, -1), 01244 ImageRef( 1, -1), 01245 01246 ImageRef(-1, 0), 01247 ImageRef( 1, 0), 01248 01249 ImageRef(-1, 1), 01250 ImageRef( 0, 1), 01251 ImageRef( 1, 1), 01252 }; 01253 for(unsigned int i=0; i < pixels.size(); i++) 01254 { 01255 if(!(log_ratios[pixels[i]] > threshold)) 01256 goto not_a_maximum; 01257 01258 for(int j=0; j < 8; j++) 01259 if(!log_ratios.in_image(pixels[i] + neighbours[j]) || ! (log_ratios[pixels[i]] > log_ratios[pixels[i] + neighbours[j]])) 01260 goto not_a_maximum; 01261 01262 spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), log_normal_mode(blur_mu, blur_sigma), pixels[i].x, pixels[i].y)); 01263 01264 //Remove the pixels around the initial spots 01265 for(unsigned int j=0; j < disc.size(); j++) 01266 valid_pixels.erase(pixels[i] + disc[j]); 01267 01268 not_a_maximum:; 01269 } 01270 01271 for(unsigned int i=0; i < spots.size(); i++) 01272 cout << spots[i] << endl; 01273 01274 //Now place down extra spots in the remaining space. 01275 while(!valid_pixels.empty()) 01276 { 01277 ImageRef p = max_element(valid_pixels.begin(), valid_pixels.end(), LessSecond())->first; 01278 spots.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), log_normal_mode(blur_mu, blur_sigma), p.x, p.y)); 01279 01280 for(unsigned int j=0; j < disc.size(); j++) 01281 valid_pixels.erase(p + disc[j]); 01282 } 01283 01284 //This line allows extra spots to be placed down around each spot already put down. 01285 //This is a shocking hack and jenerally very unpleasant. 01286 double extra_r = GV3::get<double>("extra_spots", 0, 1); 01287 vector<ImageRef> extra = getDisc(extra_r); 01288 vector<Vector<4> > more_spots; 01289 for(unsigned int i=0; i < extra.size(); i++) 01290 if(extra[i] != ImageRef_zero) 01291 for(unsigned int j=0; j < spots.size(); j++) 01292 more_spots.push_back(spots[j] + makeVector(0, 0, extra[i].x, extra[i].y) / (2*extra_r+1)); 01293 01294 copy(more_spots.begin(), more_spots.end(), back_inserter(spots)); 01295 } 01296 else 01297 { 01298 Vector<> loaded_spots = GV3::get<Vector<> >("spots.manual_spots", "", -1); 01299 01300 if(loaded_spots.size()%4 != 0) 01301 { 01302 cerr << "Loaded spot size is not a multiple of 4\n"; 01303 exit(1); 01304 } 01305 01306 else 01307 spots = spots_to_vector(loaded_spots); 01308 } 01309 01310 //Initialize the MT19937 RNG from a seed. 01311 shared_ptr<MT19937> rng(new MT19937); 01312 rng->simple_seed(GV3::get<int>("seed", 0, 1)); 01313 01314 //Load in a checkpoint (precise RNG state, iteration and pass). 01315 int start_iteration=0; 01316 int start_pass=0; 01317 if(GV3::get<bool>("checkpoint", 0, 1)) 01318 { 01319 string rng_state = GV3::get<string>("checkpoint.rng.state", "", -1); 01320 istringstream rs(rng_state); 01321 rng->read(rs); 01322 start_iteration=GV3::get<int>("checkpoint.iteration", 0, -1); 01323 start_pass=GV3::get<int>("checkpoint.pass", 0, -1); 01324 } 01325 01326 StateParameters p; 01327 p.spots = spots; 01328 p.rng = rng; 01329 p.pass = start_pass; 01330 p.iteration = start_iteration; 01331 p.pixels = pixels; 01332 01333 return p; 01334 } 01335 01336 ///Very simple and inefficient dilation function. 01337 ///@ingroup gUtility 01338 set<ImageRef> dilate_mask(const vector<ImageRef>& v, double r) 01339 { 01340 vector<ImageRef> m = getDisc(r); 01341 01342 set<ImageRef> ret; 01343 01344 for(unsigned int i=0; i < v.size(); i++) 01345 for(unsigned int j=0; j < m.size(); j++) 01346 ret.insert(v[i] + m[j]); 01347 01348 return ret; 01349 } 01350 01351 ///How far should steps in brightness be limited to? 01352 ///@ingroup gStorm 01353 double brightness_motion_limit(double mu, double sigma, bool not_one) 01354 { 01355 if(not_one) 01356 return log_normal_std(mu, sigma); 01357 else 01358 return 1; 01359 } 01360 01361 01362 01363 01364 ///Mega class which actually does the meat of the spot fitting. 01365 ///probably could be refactored a bit. 01366 ///@ingroup gStorm 01367 class FitSpots 01368 { 01369 const vector<Image<float> >& ims; ///< Input data 01370 FitSpotsGraphics& graphics; ///< Graphics class. 01371 UserInterfaceCallback& ui; ///< Callbacks to provide user interface 01372 const vector<ImageRef> pixels; ///< Area in which to perform model fitting 01373 01374 //Spot positions 01375 vector<Vector<4> > spots; ///< State in terms of current spot positions 01376 01377 //Starting point 01378 const int start_iteration; ///< Starting iteration number (for restarting from checkpoint) 01379 int start_pass; ///< Starting pass (for restarting from checkpoint) 01380 01381 MT19937& rng; ///< Random numbewr generator 01382 01383 const double variance; ///< Variance of noise in data. Must be 1. 01384 const double intensity_mu; ///< Prior for spot intensity 01385 const double intensity_sigma; ///< Prior for spot intensity 01386 const double blur_mu; ///< Prior for spot shape 01387 const double blur_sigma; ///< Prior for spt shape 01388 01389 //pixels is dilated slightly, by a disk of size area_extra_radius. 01390 //The result is stored in allowed_area. This is used to implement a uniform 01391 //prior over position, which is uniform within allowed_area, and zero 01392 //everywhere else. 01393 // 01394 //This is implemented as a set, because the mask may extend beyond the borders of the 01395 //image slightly 01396 const double area_extra_radius; ///< Extra size beyone marked region in which spots are allowed to exist 01397 set<ImageRef> allowed_area; ///< Total allowed region, pixels dilated by area_extra_radius 01398 const int use_position_prior; ///< Should a proper prior over position be uesd? A clue: yes. 01399 const double position_prior; ///< Value for the posision prior, i.e. reciprocal of area 01400 01401 //General optimizing and sampling parameters 01402 const double max_motion; ///< Maximum motion on any axes for optimization. See ConjugateGradientOnly. 01403 const int sample_iterations; ///< Number of mixing samples to use in Gibbs sampling 01404 01405 //Task specific optimizing and sampling parameters: 01406 01407 //Main optimization loop 01408 const int main_cg_max_iterations; ///< Maximum iterations allowed for ConjugateGradientOnly in main optimization loop 01409 const int main_samples; ///< Number of samples to use in main loop 01410 const int main_passes; ///< Number of passes to perform per iteration in main loop 01411 const int outer_loop_iterations; ///< Total number of iterations to perform 01412 01413 //Spot selection loop 01414 const int add_remove_tries; ///< Number of attemts to add/remove a spot 01415 const int add_remove_opt_samples; ///< Number of samples to use in model modification phase 01416 const int add_remove_opt_retries; ///< Number of attempts restarting the optimization to escape saddle points 01417 const int add_remove_opt_hess_inner_samples; ///< Number of extra FFBS samples to use for computing Hessian when testing for convergence to non-saddle point 01418 const int h_outer_samples; ///< Number of samples used for computing Hessian as part of Laplace's approximation 01419 const int h_inner_samples; ///< Number of additional FFBS samples to use for computing Hessian as part of Laplace's approximation 01420 const int tsamples; ///< Number of samples to use in thermodynamic integration 01421 01422 // Average of all inputs (used for debugging) 01423 const Image<float> ave; ///< Average of input data: used for 01424 01425 01426 //File to save output to 01427 ofstream& save_spots; ///< Output stream for log file 01428 01429 //Time accumulators for benchmarking 01430 double time_gibbs; ///< Benchmarking data 01431 double time_cg; ///< Benchmarking data 01432 01433 ///Motion limit for ConjugateGradientOnly 01434 ///The x distance, y distance and spot size are all approximately the same scale 01435 ///which is of order 1. The brigntness is not. By default, the brightness limit is also 1. 01436 ///This flag controls whether is should be set to the standard deviation of the brightness prior 01437 ///distribution. This setting will put the motion limit on the same scale as the other three 01438 ///parameters. 01439 const bool scale_brightness_limit; 01440 const Vector<4> limit; ///< Limit vector for ConjugateGradientOnly 01441 01442 const Matrix<3> A; ///< Transition matrix 01443 const Vector<3> pi; ///< Initial probabilities 01444 01445 01446 vector<vector<double> > pixel_intensities; ///<Pixel intensities for all images [frame][pixel] 01447 01448 DataForMCMC data_for_t_mcmc; ///<Aggergated data for thermodynamic integration 01449 DataForMCMC data_for_h_mcmc; ///<Aggergated data for finding hessian 01450 01451 int iteration; ///< Iteration number 01452 01453 01454 TIME(cvd_timer timer;) 01455 public: 01456 01457 FitSpots(const vector<Image<float> >& ims_, FitSpotsGraphics& graphics_, UserInterfaceCallback& ui_, StateParameters& params, ofstream& save_spots_) 01458 :ims(ims_), graphics(graphics_),ui(ui_), 01459 01460 //Set up the main parameters 01461 pixels(params.pixels), 01462 spots(params.spots), 01463 start_iteration(params.iteration), 01464 start_pass(params.pass), 01465 rng(*params.rng), 01466 01467 01468 01469 //Set up all the system parameters 01470 variance(1), // it should be 01471 01472 //To scale the X axis of a log-normal distribution, only 01473 //the mu parameter needs to be changed... 01474 intensity_mu(GV3::get<double>("intensity.rel_mu", 0., -1) + log(sqrt(variance))), 01475 intensity_sigma(GV3::get<double>("intensity.rel_sigma", 0., -1)), 01476 blur_mu(GV3::get<double>("blur.mu", 0., -1)), 01477 blur_sigma(GV3::get<double>("blur.sigma", 0., -1)), 01478 01479 //Spot position prior 01480 area_extra_radius(GV3::get<double>("position.extra_radius", 0., -1)), 01481 allowed_area(dilate_mask(pixels, area_extra_radius)), 01482 use_position_prior(GV3::get<bool>("position.use_prior", true, -1)), 01483 position_prior(1.0 / allowed_area.size()), 01484 01485 //General optimizing and sampling parameters 01486 max_motion(GV3::get<double>("cg.max_motion", 0., -1)), 01487 sample_iterations(GV3::get<int>("gibbs.mixing_iterations", 0, -1)), 01488 01489 //Task specific optimizing and sampling parameters: 01490 01491 //Main optimization loop 01492 main_cg_max_iterations(GV3::get<double>("main.cg.max_iterations", 0., -1)), 01493 main_samples(GV3::get<int>("main.gibbs.samples", 0, -1)), 01494 main_passes(GV3::get<int>("main.passes", 0, -1)), 01495 outer_loop_iterations(GV3::get<int>("main.total_iterations", 100000000, 1)), 01496 01497 //Spot selection loop 01498 add_remove_tries(GV3::get<int>("add_remove.tries", 0, -1)), 01499 add_remove_opt_samples(GV3::get<int>("add_remove.optimizer.samples", 0, -1)), 01500 add_remove_opt_retries(GV3::get<int>("add_remove.optimizer.attempts", 0, -1)), 01501 add_remove_opt_hess_inner_samples(GV3::get<int>("add_remove.optimizer.hessian_inner_samples", 0, -1)), 01502 h_outer_samples(GV3::get<int>("add_remove.hessian.outer_samples", 0, -1)), 01503 h_inner_samples(GV3::get<int>("add_remove.hessian.inner_samples", 0, -1)), 01504 tsamples(GV3::get<int>("add_remove.thermo.samples", 0, -1)), 01505 01506 ave(average_image(ims_)), 01507 01508 save_spots(save_spots_), 01509 01510 time_gibbs(0), 01511 time_cg(0), 01512 01513 scale_brightness_limit(GV3::get<bool>("max_motion.use_brightness_std", 0, -1)), 01514 limit(makeVector(brightness_motion_limit(intensity_mu, intensity_sigma, scale_brightness_limit), 1, 1, 1)*max_motion), 01515 01516 A(GV3::get<Matrix<3> >("A", Zeros, 1)), 01517 pi(GV3::get<Vector<3> >("pi", Zeros, 1)), 01518 01519 01520 data_for_t_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, tsamples, sample_iterations, A, pi, rng), 01521 data_for_h_mcmc(pixels, pixel_intensities, intensity_mu, intensity_sigma, blur_mu, blur_sigma, variance, h_outer_samples, sample_iterations, A, pi, rng) 01522 01523 { 01524 assert_same_size(ims); 01525 01526 //Pixel intensities for all images [frame][pixel] 01527 pixel_intensities.resize(ims.size(), vector<double>(pixels.size())); 01528 for(unsigned int frame=0; frame < ims.size(); frame++) 01529 for(unsigned int p=0; p < pixels.size(); p++) 01530 pixel_intensities[frame][p] = ims[frame][pixels[p]]; 01531 01532 } 01533 01534 ///Perform a complete iteration of the optimzation stage of the spot firrint algorithm 01535 void optimize_each_spot_in_turn_for_several_passes() 01536 { 01537 //Precompute the intensities for all spot pixels 01538 vector<vector<double> > spot_intensities; //[spot][pixel] 01539 for(unsigned int i=0; i < spots.size(); i++) 01540 spot_intensities.push_back(compute_spot_intensity(pixels, spots[i])); 01541 01542 //Which pixels does each spot have? 01543 vector<vector<int> > spot_pixels; //[spot][pixel] 01544 spot_pixels.resize(spots.size()); 01545 for(unsigned int s=0; s < spots.size(); s++) 01546 get_spot_pixels(pixels, spots[s], spot_pixels[s]); 01547 01548 01549 //Optimize the model, N spots at a time. 01550 // 01551 for(int pass=start_pass; pass < main_passes; pass++) 01552 { 01553 save_spots << "Pass: " << pass << endl; 01554 rng.write(save_spots); 01555 save_spots << endl; 01556 01557 start_pass=0; // This is only nonzero first time, since we might chekpoint midway through an iteration 01558 save_spots << "PASS" << pass << ": " << setprecision(20) << scientific << spots_to_Vector(spots) << endl; 01559 save_spots << "ENDCHECKPOINT" << endl << flush; 01560 01561 ui.per_pass(iteration, pass, spots); 01562 //Sort spots according to pass%4 01563 01564 //Sort the spots so that the optimization runs in sweeps 01565 //This heiristic seems to increase the rate of propagation of information 01566 //about spot positions. 01567 01568 //Create a list of indices 01569 vector<int> index = sequence(spots.size()); 01570 01571 int passs = pass + iteration; 01572 //Sort the indices according to the position of the spot that they index 01573 if(passs%4 == 0) 01574 sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 2>(spots)); 01575 else if(passs%4==1) 01576 sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 2>(spots)); 01577 else if(passs%4==1) 01578 sort(index.begin(), index.end(), IndexLexicographicPosition<less<double>, 3>(spots)); 01579 else 01580 sort(index.begin(), index.end(), IndexLexicographicPosition<greater<double>, 3>(spots)); 01581 01582 //Reorder the spots and their intensities and their pixel lists 01583 { 01584 vector<Vector<4> > tmp_spot(index.size()); 01585 vector<vector<double> > tmp_spot_intensities(index.size()); 01586 vector<vector<int> > tmp_spot_pixels(index.size()); 01587 for(unsigned int i=0; i < index.size(); i++) 01588 { 01589 tmp_spot[i] = spots[index[i]]; 01590 swap(tmp_spot_intensities[i], spot_intensities[index[i]]); 01591 swap(tmp_spot_pixels[i], spot_pixels[i]); 01592 } 01593 01594 swap(tmp_spot, spots); 01595 swap(tmp_spot_intensities, spot_intensities); 01596 swap(tmp_spot_pixels, spot_pixels); 01597 } 01598 01599 //Sweep through and optimize each spot in turn 01600 for(int s=0; s < (int)spots.size(); s++) 01601 { 01602 ui.per_spot(iteration, pass, s, spots.size()); 01603 ui.perhaps_stop(); 01604 01605 TIME(timer.reset();) 01606 //Get some samples with Gibbs sampling 01607 vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling 01608 vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel] 01609 01610 GibbsSampler2 sampler(pixel_intensities, spot_intensities, spots, spot_pixels, A, pi, variance, sample_iterations); 01611 for(int i=0; i < main_samples; i++) 01612 { 01613 sampler.next(rng); 01614 sample_list.push_back(sampler.sample()); 01615 sample_intensities.push_back(sampler.sample_intensities()); 01616 01617 ui.perhaps_stop(); 01618 } 01619 01620 //First, remove the spot from all the samples. 01621 for(unsigned int i=0; i < sample_list.size(); i++) 01622 remove_spot(sample_intensities[i], spot_intensities[s], sample_list[i][s]); 01623 01624 //cout << timer.get_time() << endl; 01625 TIME(time_gibbs += timer.reset();) 01626 01627 //Package up all the data 01628 SampledBackgroundData data(sample_intensities, pixel_intensities, pixels, 01629 intensity_mu, intensity_sigma, blur_mu, blur_sigma, 01630 A, pi, variance); 01631 01632 //Derivative computer: 01633 SpotNegProbabilityDiffWithSampledBackground compute_deriv(data); 01634 01635 01636 graphics.draw_pixels(pixels, 0, 0, 1, 1); 01637 graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s); 01638 graphics.swap(); 01639 01640 //Optimize spot "s" 01641 //Optimize with the derivatives only since the actual probability 01642 //is much harder to compute 01643 ConjugateGradientOnly<4> cg(spots[s], compute_deriv, limit); 01644 01645 01646 cg.max_iterations = main_cg_max_iterations; 01647 01648 01649 #if 0 01650 cout << setprecision(10); 01651 cout << spots_to_Vector(spots) << endl; 01652 Matrix<4> hess, hess_errors; 01653 cout << "Hello, my name is Inigo Montoya\n"; 01654 /*for(int i=0; i < 4; i++) 01655 { 01656 Matrix<4, 2> d = numerical_gradient_with_errors(NthDeriv(compute_deriv, i), cg.x); 01657 hess[i] = d.T()[0]; 01658 hess_errors[i] = d.T()[1]; 01659 } 01660 */ 01661 //cout << "Errors:\n" << hess_errors << endl; 01662 //cout << "NHess:\n" << hess<< endl; 01663 Matrix<4> rhess = -sampled_background_spot_hessian(cg.x, data); 01664 cout << "Hess:\n" << rhess << endl; 01665 //cout << "Err:\n" << hess - rhess << endl; 01666 01667 //Vector<4> grad = compute_deriv(cg.x); 01668 01669 //Matrix<4> e = hess - rhess; 01670 01671 //for(int i=0; i < 4; i++) 01672 // for(int j=0; j < 4; j++) 01673 // e[i][j] /= hess_errors[i][j]; 01674 01675 //cout << "Err:\n" << e << endl; 01676 cout << "Deriv:" << compute_deriv(cg.x) << endl; 01677 //cout << "Full:\n" << sampled_background_spot_hessian2(cg.x, data) - grad.as_col()*grad.as_row() << endl; 01678 01679 FreeEnergyHessian hesscomputer(data_for_h_mcmc); 01680 01681 Matrix<4> nhess = hesscomputer.hessian(spots, 0); 01682 cout << "NHess:\n" << nhess << endl; 01683 01684 cout << "Turbo-N Hess:\n" << sampled_background_spot_hessian_ffbs(cg.x, data, 10000) << endl; 01685 01686 cout << "TI energy: " << NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(spots)) << endl; 01687 cout << "FA energy: " << SpotProbabilityWithSampledBackgroundFAKE(data)(cg.x) << endl; 01688 01689 01690 01691 //cout << "Numerical hessian from FD:\n" << numerical_hessian(SpotProbabilityWithSampledBackgroundFAKE(data), cg.x) << endl; 01692 exit(0); 01693 #endif 01694 //cout << "Starting opt... " << cg.x << endl; 01695 while(cg.iterate(compute_deriv)) 01696 { 01697 graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), s, cg.x); 01698 graphics.draw_pixels(pixels, 0, 0, 1, .2); 01699 graphics.swap(); 01700 ui.perhaps_stop(); 01701 //cout << cg.x << endl; 01702 } 01703 01704 //Update to use the result of the optimization 01705 spots[s] = cg.x; 01706 //cout << "End: " << cg.x << endl; 01707 01708 graphics.draw_krap(spots, scale_to_bytes(ave), boundingbox(pixels), -1); 01709 graphics.swap(); 01710 01711 //Recompute the new spot intensity, since the spot has changed 01712 spot_intensities[s] = compute_spot_intensity(pixels, spots[s]); 01713 01714 //Recompute which are the useful pixels 01715 get_spot_pixels(pixels, spots[s], spot_pixels[s]); 01716 01717 //Is the spot within the allowed area, i.e. is it's prior 0? 01718 //The prior is sero only if it we are using it and we're in an invalid area 01719 01720 ImageRef quantized_spot_position = ir_rounded(spots[s].slice<2,2>()); 01721 bool zero_prior = use_position_prior && (allowed_area.count(quantized_spot_position)==0); 01722 01723 //Check to see if spot has been ejected. If spot_pixels is empty, then it has certainly been ejected. 01724 if(spot_pixels[s].empty() || zero_prior) 01725 { 01726 //Spot has been ejected. Erase it. 01727 cout << " Erasing ejected spot: " << spot_pixels[s].empty() << " " << zero_prior << endl; 01728 cout << spots[s] << endl; 01729 //GUI_Pause(1); 01730 01731 spot_intensities.erase(spot_intensities.begin() + s); 01732 spot_pixels.erase(spot_pixels.begin() + s); 01733 spots.erase(spots.begin() + s); 01734 s--; 01735 //exit(0); 01736 } 01737 01738 //cout << timer.get_time() << endl; 01739 TIME(time_cg += timer.reset();) 01740 01741 //cout << "Times: " << time_gibbs << " " << time_cg << endl; 01742 //save_spots << "INTERMEDIATE: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl; 01743 } 01744 } 01745 } 01746 01747 ///Perform a complete iteration of the model size modification stage of the spot fitting algorithm 01748 void try_modifying_model() 01749 { 01750 01751 for(int i=0; i < add_remove_tries; i++) 01752 { 01753 ui.per_modification(iteration, i, add_remove_tries); 01754 ui.perhaps_stop(); 01755 01756 cout << endl << endl << "Modifying the model" << endl << "======================\n"; 01757 cout << "Hello\n"; 01758 bool add_spot = (rng() > 0.5) || (spots.size() == 1); 01759 cout << "World\n"; 01760 01761 vector<Vector<4> > model_1, model_2; 01762 01763 01764 int original; //What is the original model? Model 1 or Model 2? 01765 01766 if(add_spot) 01767 { 01768 model_1 = spots; 01769 model_2 = model_1; 01770 01771 //Pick a pixel within the thresholded ones as a starting point 01772 int r; 01773 do 01774 { 01775 r = (int)(rng() * pixels.size()); 01776 //cout << r << " " << log_ratios[pixels[r]] << " " << pixels[r] << " " << threshold << endl; 01777 } 01778 while(0); 01779 01780 //This version does not (yet?) suppotrt thresholding on log_ratios 01781 //for placing spots, since the purpose of log_ratios has diminished. 01782 //while(log_ratios[pixels[r]] < threshold); 01783 01784 01785 model_2.push_back(makeVector(log_normal_mode(intensity_mu, intensity_sigma), 01786 log_normal_mode(blur_mu, blur_sigma), 01787 pixels[r].x + rng()-.5, pixels[r].y + rng() - .5)); 01788 cout << "Adding a spot\n"; 01789 01790 original = 1; 01791 } 01792 else 01793 { 01794 //Pick a random point to optimize/remove 01795 int a_random_spot = static_cast<int>(rng() * spots.size()); 01796 model_1 = spots; 01797 swap(model_1[model_1.size()-1], model_1[a_random_spot]); 01798 01799 model_2 = model_1; 01800 01801 model_1.pop_back(); 01802 cout << "Removing a spot\n"; 01803 original = 2; 01804 } 01805 01806 //The mobile spot is always the last spot of model_2 01807 const int spot = model_2.size() - 1; 01808 01809 cout << "Original model: " << original << endl; 01810 01811 //Precompute the intensities for all spot pixels 01812 //Model 2 is always one longer than model 1 and 01813 //differs only on the extra element 01814 vector<vector<double> > model2_spot_intensities; //[spot][pixel] 01815 for(unsigned int i=0; i < model_2.size(); i++) 01816 model2_spot_intensities.push_back(compute_spot_intensity(pixels, model_2[i])); 01817 01818 //Which pixels does each spot have? 01819 vector<vector<int> > model2_spot_pixels(model_2.size()); //[spot][pixel] 01820 for(unsigned int s=0; s < model_2.size(); s++) 01821 get_spot_pixels(pixels, model_2[s], model2_spot_pixels[s]); 01822 01823 //Optimize spot: 01824 { 01825 cout << "Optimizing spot for model selection\n"; 01826 01827 01828 //Get some samples with Gibbs sampling 01829 vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling 01830 vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel] 01831 01832 GibbsSampler2 sampler(pixel_intensities, model2_spot_intensities, model_2, model2_spot_pixels, A, pi, variance, sample_iterations); 01833 for(int i=0; i < add_remove_opt_samples; i++) 01834 { 01835 sampler.next(rng); 01836 sample_list.push_back(sampler.sample()); 01837 sample_intensities.push_back(sampler.sample_intensities()); 01838 ui.perhaps_stop(); 01839 } 01840 01841 //First, remove the spot from all the samples. 01842 for(unsigned int i=0; i < sample_list.size(); i++) 01843 remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]); 01844 01845 //Package up all the data 01846 SampledBackgroundData data(sample_intensities, pixel_intensities, pixels, 01847 intensity_mu, intensity_sigma, blur_mu, blur_sigma, 01848 A, pi, variance); 01849 01850 //Derivative computer: 01851 SpotNegProbabilityDiffWithSampledBackground compute_deriv(data); 01852 01853 graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot); 01854 graphics.swap(); 01855 01856 //Optimize spot "s" 01857 //Optimize with the derivatives only since the actual probability 01858 //is much harder to compute 01859 ConjugateGradientOnly<4> cg(model_2[spot], compute_deriv, limit); 01860 01861 01862 for(int attempt=0; attempt < add_remove_opt_retries; attempt++) 01863 { 01864 cout << "Attempt " << attempt << " of " << add_remove_opt_retries << endl; 01865 ui.perhaps_stop(); 01866 01867 //Optimize with conjugate gradient 01868 while(cg.iterate(compute_deriv)) 01869 { 01870 ui.perhaps_stop(); 01871 graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), spot, cg.x); 01872 graphics.swap(); 01873 } 01874 01875 //Check for being at a saddle point (no point checking on the last try) 01876 //All eigenvectors should be negative at a maximum. 01877 //WTF: is this a bug? WTF? 01878 //It was this: 01879 //if(attempt < add_remove_tries - 1) 01880 if(attempt < add_remove_opt_retries - 1) 01881 { 01882 Matrix<4> hessian = sampled_background_spot_hessian_ffbs(cg.x, data, add_remove_opt_hess_inner_samples, rng); 01883 SymEigen<4> hess_decomp(hessian); 01884 01885 //cout << "What the fuck:" << endl << hessian << endl << hessian3<< endl << hessian2 << endl; 01886 01887 cout << "Eigenvalues are: " << hess_decomp.get_evalues() << endl; 01888 01889 if(hess_decomp.is_negdef()) 01890 break; 01891 else 01892 { 01893 //Restart in the direction of the best uphill part 01894 cg.init(cg.x + 0.1 * hess_decomp.get_evectors()[3], (hess_decomp.get_evectors()[3])); 01895 01896 01897 cout << "Grad = " << compute_deriv(cg.x) << endl; 01898 for(int i=0; i < 4; i++) 01899 { 01900 cout << "Direction: " << i << endl; 01901 cout << unit(compute_deriv(cg.x + 0.1*hess_decomp.get_evectors()[i])) * hess_decomp.get_evectors()[i] << endl; 01902 } 01903 01904 for(int i=0; i < 4; i++) 01905 { 01906 cout << "Direction: " << i << endl; 01907 Vector<4> d = Zeros; 01908 d[i] = 1; 01909 cout << unit(compute_deriv(cg.x + d)) * hess_decomp.get_evectors()[i] << endl; 01910 } 01911 } 01912 } 01913 } 01914 01915 01916 //Update to use the result of the optimization 01917 model_2[spot] = cg.x; 01918 01919 graphics.draw_krap(model_2, scale_to_bytes(ave), boundingbox(pixels), -1); 01920 graphics.swap(); 01921 01922 01923 //Update cached data based on the new spot position 01924 model2_spot_intensities[spot] = compute_spot_intensity(pixels, model_2[spot]); 01925 get_spot_pixels(pixels, model_2[spot], model2_spot_pixels[spot]); 01926 01927 cout << "Done optimizing for model selection\n"; 01928 } 01929 01930 01931 //Which model to keep? 01932 int keep=original; 01933 01934 //Compute position prior (and we might be able to reject it really quickly here) 01935 bool zero_prior = use_position_prior && (allowed_area.count(ir_rounded(model_2[spot].slice<2,2>()))==0); 01936 01937 if(zero_prior) 01938 { 01939 //Model 2 went bad, so we clearly keep model 1 01940 keep = 1; 01941 } 01942 else 01943 { 01944 01945 //The position prior is independent 01946 //Compute the difference model2 - model1 01947 //This is only valid if model2 is in the valid region 01948 double position_log_prior_model2_minus_model1; 01949 if(use_position_prior) 01950 position_log_prior_model2_minus_model1 = (model_2.size() - model_1.size()) * ln(position_prior); 01951 else 01952 position_log_prior_model2_minus_model1 = 0; 01953 01954 01955 //Integrate model_2 01956 //First compute the Hessian since this might go wrong. 01957 01958 //FreeEnergyHessian hesscomputer(data_for_h_mcmc); 01959 Matrix<4> hess;// = hesscomputer.hessian(model_2, spot); 01960 01961 //Use turbohess here since it is much faster, as the backwards sampling step is fast 01962 //We expect this hessian to be really quite different, if the spot has moved from 01963 //a long way away, since the sampling will have changed dramatically 01964 { 01965 //Get some samples with Gibbs sampling 01966 vector<vector<vector<State> > > sample_list; //[N][spot][frame]: list of samples drawn using Gibbs sampling 01967 vector<vector<vector<double> > > sample_intensities; //[sample][frame][pixel] 01968 01969 GibbsSampler sampler(pixel_intensities, model2_spot_intensities, model_2, A, pi, variance, sample_iterations); 01970 for(int i=0; i < h_outer_samples; i++) 01971 { 01972 ui.perhaps_stop(); 01973 sampler.next(rng); 01974 sample_list.push_back(sampler.sample()); 01975 sample_intensities.push_back(sampler.sample_intensities()); 01976 } 01977 01978 //First, remove the spot from all the samples. 01979 for(unsigned int i=0; i < sample_list.size(); i++) 01980 remove_spot(sample_intensities[i], model2_spot_intensities[spot], sample_list[i][spot]); 01981 01982 //Package up all the data 01983 SampledBackgroundData data(sample_intensities, pixel_intensities, pixels, 01984 intensity_mu, intensity_sigma, blur_mu, blur_sigma, 01985 A, pi, variance); 01986 01987 hess = sampled_background_spot_hessian_ffbs(model_2[spot], data, h_inner_samples, rng); 01988 } 01989 01990 01991 double det = determinant(-hess / (M_PI*2)); 01992 SymEigen<4> hess_decomp(-hess); 01993 cout << "Hessien Eigenvalues are: " << hess_decomp.get_evalues() << endl; 01994 const double smallest_evalue = 1e-6; 01995 01996 //If the determinant is negative, then we are still at a saddle point 01997 //despite the attempts above, so abandon the attempt. 01998 if(hess_decomp.get_evalues()[0] > smallest_evalue) 01999 { 02000 02001 //Compute the free energy of model 2 at the MLE estimate 02002 cout << "Model 2:\n"; 02003 // double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_2)); 02004 double model_2_energy = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_2), model2_spot_pixels); 02005 cout << "Energy: " << model_2_energy << endl; 02006 02007 //Combine the MLE energy and Hessian using Laplace's approximation 02008 double model_2_occam = -0.5 * log(det); 02009 double model_2_prob = model_2_energy + model_2_occam + position_log_prior_model2_minus_model1; 02010 02011 cout << "Occam: " << model_2_occam << endl; 02012 cout << "Position: " << position_log_prior_model2_minus_model1 << endl; 02013 cout << "Prob: " << model_2_prob << endl; 02014 02015 //Integrate model_1 02016 //It has no parameters, in this formulation. 02017 //double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1)); 02018 02019 //Note that model_1 always has one fewer spots, and the last spot is always the one 02020 //missing, so we can make the correct mask very easily: 02021 model2_spot_pixels.pop_back(); 02022 double model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc).compute_with_mask(spots_to_Vector(model_1), model2_spot_pixels); 02023 cout << "Prob: " << model_1_prob << endl; 02024 //model_1_prob = -NegativeFreeEnergy(data_for_t_mcmc)(spots_to_Vector(model_1)); 02025 02026 if(model_2_prob > model_1_prob) 02027 keep=2; 02028 else 02029 keep=1; 02030 02031 cout << "Models evaluated\n"; 02032 } 02033 else 02034 { 02035 cout << "Determinant has bad eigenvalues!\n"; 02036 keep = original; 02037 cout << hess_decomp.get_evalues() << endl; 02038 } 02039 } 02040 02041 if(keep == 2) 02042 { 02043 spots = model_2; 02044 cout << "Keeping model 2\n"; 02045 02046 } 02047 else 02048 { 02049 spots = model_1; 02050 cout << "Keeping model 1\n"; 02051 } 02052 02053 if(original != keep) 02054 { 02055 cout << "Model changed!\n"; 02056 //break; 02057 } 02058 } 02059 } 02060 02061 ///Run the complete optimization algorithm. 02062 void run() 02063 { 02064 graphics.init(ims[0].size()); 02065 save_spots << "LOGVERSION 2" << endl; 02066 save_spots << "PIXELS"; 02067 for(unsigned int i=0; i < pixels.size(); i++) 02068 save_spots << " " << pixels[i].x << " " << pixels[i].y; 02069 save_spots << endl; 02070 02071 save_spots << "BEGINGVARLIST" << endl; 02072 GV3::print_var_list(save_spots, "", 1); 02073 save_spots << "ENDGVARLIST" << endl; 02074 02075 //TODO all GVARS are set, so dump out gvars. 02076 02077 cout << "Limit vector: " << limit << endl; 02078 02079 for(iteration=start_iteration; iteration < outer_loop_iterations ;iteration++) 02080 { 02081 save_spots << "Iteration: " << iteration << " (" << iteration * main_passes << ")" << endl; 02082 save_spots << "MAIN: " << setprecision(20) << scientific << spots_to_Vector(spots) << endl; 02083 02084 cout << endl << endl << "----------------------" << endl << "Optimizing:\n"; 02085 cout << spots.size() << endl; 02086 02087 02088 optimize_each_spot_in_turn_for_several_passes(); 02089 02090 //spot_intensities is be correct here! 02091 try_modifying_model(); 02092 } 02093 save_spots << "FINAL: " << setprecision(15) << scientific << spots_to_Vector(spots) << endl; 02094 } 02095 02096 }; 02097 02098 ///Wrapper function for using FitSpots 02099 ///@ingroup gStorm 02100 void fit_spots_new(const vector<Image<float> >& ims, StateParameters& p, ofstream& save_spots, FitSpotsGraphics& gr) 02101 { 02102 auto_ptr<UserInterfaceCallback> ui = null_ui(); 02103 FitSpots fit(ims, gr, *ui, p, save_spots); 02104 fit.run(); 02105 } 02106 02107 ///Wrapper function for using FitSpots 02108 ///@ingroup gStorm 02109 void fit_spots_new(const vector<Image<float> >& ims, StateParameters& p, ofstream& save_spots, FitSpotsGraphics& gr, UserInterfaceCallback& ui) 02110 { 02111 try{ 02112 FitSpots fit(ims, gr, ui, p, save_spots); 02113 fit.run(); 02114 } 02115 catch(UserInterfaceCallback::UserIssuedStop) 02116 { 02117 } 02118 }