-
Notifications
You must be signed in to change notification settings - Fork 27
/
simulate.cpp
3013 lines (2476 loc) · 157 KB
/
simulate.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
#include "../../core/core_headers.h"
#include "../../core/water.h"
#include "../../constants/constants.h"
#include "calc_hydration_radius.h"
#include "../../core/scattering_potential.h"
#include "wave_function_propagator.h"
#include "../../core/padded_coordinates.h"
const int MAX_3D_SIZE = 1536; // ~14.5 Gb in single precision
using namespace cistem_timer;
//#define DO_BEAM_TILT_FULL true
const float DISTANCE_INIT = 100000.0f; // Set the distance slab to a large value
// I would have thought this should be 1/4, however, when calculating a 3D and then plotting ln(F) vs freq^-2, it is clear that the relative difference
// requires this factor to be 1/2. This is somewhat bothersome and ideally would have an explanation.
const float BfactorFactor = 0.5f;
const int N_WATER_TERMS = 40;
const float WATER_BFACTOR_PER_ELECTRON_PER_SQANG = 34.0f;
const float PHASE_PLATE_BFACTOR = 4.0f;
const float inelastic_scalar_water = 0.0725f; // this value is for 1 Apix. When not taking the sqrt (original but I think wrong approach) the value 0.75 worked
// Parameters for calculating water. Because all waters are the same, except a sub-pixel origin offset. A SUB_PIXEL_NeL stack of projected potentials is pre-calculated with given offsets.
const AtomType SOLVENT_TYPE = oxygen; // 2 N, 3 O 15, fit water
const float water_oxygen_ratio = 1.0f; //0.8775;
const int SUB_PIXEL_NEIGHBORHOOD = 2;
const int SUB_PIXEL_NeL = pow((SUB_PIXEL_NEIGHBORHOOD * 2 + 1), 3);
const float TAPER[29] = {0, 0, 0, 0, 0,
0.003943, 0.015708, 0.035112, 0.061847, 0.095492, 0.135516, 0.181288, 0.232087,
0.287110, 0.345492, 0.406309, 0.468605, 0.531395, 0.593691, 0.654508, 0.712890,
0.767913, 0.818712, 0.864484, 0.904508, 0.938153, 0.964888, 0.984292, 0.996057};
const int IMAGEPADDING = 0; //512; // Padding applied for the Fourier operations in the propagation steps. Removed prior to writing out.
const int IMAGETRIMVAL = 0; //IMAGEPADDING + 2*TAPERWIDTH;
const float MAX_PIXEL_SIZE = 4.0f;
//beam_tilt_xform
//const int n_tilt_angles = 101;
//const float SET_TILT_ANGLES[n_tilt_angles] = {-70.000, -68.600, -67.200, -65.800, -64.400, -63.000, -61.600, -60.200, -58.800, -57.400, -56.000, -54.600, -53.200, -51.800, -50.400, -49.000, -47.600, -46.200, -44.800, -43.400, -42.000, -40.600, -39.200, -37.800, -36.400, -35.000, -33.600, -32.200, -30.800, -29.400, -28.000, -26.600, -25.200, -23.800, -22.400, -21.000, -19.600, -18.200, -16.800, -15.400, -14.000, -12.600, -11.200, -9.800, -8.400, -7.000, -5.600, -4.200, -2.800, -1.400, 0.000, 1.400, 2.800, 4.200, 5.600, 7.000, 8.400, 9.800, 11.200, 12.600, 14.000, 15.400, 16.800, 18.200, 19.600, 21.000, 22.400, 23.800, 25.200, 26.600, 28.000, 29.400, 30.800, 32.200, 33.600, 35.000, 36.400, 37.800, 39.200, 40.600, 42.000, 43.400, 44.800, 46.200, 47.600, 49.000, 50.400, 51.800, 53.200, 54.600, 56.000, 57.400, 58.800, 60.200, 61.600, 63.000, 64.400, 65.800, 67.200, 68.600, 70.000};
//const float SET_TILT_ANGLES[n_tilt_angles] = {-15.400, -14.000, -12.600, -11.200, -9.800, -8.400, -7.000, -5.600, -4.200, -2.800, -1.400, -0.000, 1.400, 2.800, 4.200, 5.600, 7.000, 8.400, 9.800, 11.200, 12.600, 14.000, 15.400, 16.800, 18.200, 19.600, 21.000, 22.400, 23.800, 25.200, 26.600, 28.000, 29.400, 30.800, 32.200, 33.600, 35.000, 36.400, 37.800, 39.200, 40.600, 42.000, 43.400, 44.800, 46.200, 47.600, 49.000, 50.400, 51.800, 53.200, 54.600, 56.000, 57.400, 58.800, 60.200, 61.600, 63.000, 64.400, 65.800, 67.200, 68.600, 70.000, -16.800, -18.200, -19.600, -21.000, -22.400, -23.800, -25.200, -26.600, -28.000, -29.400, -30.800, -32.200, -33.600, -35.000, -36.400, -37.800, -39.200, -40.600, -42.000, -43.400, -44.800, -46.200, -47.600, -49.000, -50.400, -51.800, -53.200, -54.600, -56.000, -57.400, -58.800, -60.200, -61.600, -63.000, -64.400, -65.800, -67.200, -68.600, -70.000};
//const int n_tilt_angK2Count-300kv-50e-3eps.txtles = 50;
//const float SET_TILT_ANGLES[n_tilt_angles] = {-15.400, -12.600, -9.800, -7.000, -4.200, -1.400, 1.400, 4.200, 7.000, 9.800, 12.600, 15.400, 18.200, 21.000, 23.800, 26.600, 29.400, 32.200, 35.000, 37.800, 40.600, 43.400, 46.200, 49.000, 51.800, 54.600, 57.400, 60.200, 63.000, 65.800, 68.600, -18.200, -21.000, -23.800, -26.600, -29.400, -32.200, -35.000, -37.800, -40.600, -43.400, -46.200, -49.000, -51.800, -54.600, -57.400, -60.200, -63.000, -65.800, -68.600};
//const int n_tilt_angles = 51;
//const float SET_TILT_ANGLES[n_tilt_angles] = {0.00, 2.800, -2.800, 5.600, -5.600, 8.400, -8.400, 11.200, -11.200, 14.000, -14.000, 16.800, -16.800, 19.600, -19.600, 22.400, -22.400, 25.200, -25.200, 28.000, -28.000, 30.800, -30.800, 33.600, -33.600, 36.400, -36.400, 39.200, -39.200, 42.000, -42.000, 44.800, -44.800, 47.600, -47.600, 50.400, -50.400, 53.200, -53.200, 56.000, -56.000, 58.800, -58.800, 61.600, -61.600, 64.400, -64.400, 67.200, -67.200, 70.000, -70.000};
// Some of the more common elements, should add to this later. These are from Peng et al. 1996.
//const int n_tilt_angles = 41;
//const float SET_TILT_ANGLES[n_tilt_angles] = {0, 3, -3, 6, -6, 9, -9, 12, -12, 15, -15, 18, -18, 21, -21, 24, -24, 27, -27, 30, -30, 33, -33, 36, -36, 39, -39, 42, -42, 45, -45, 48, -48, 51, -51, 54, -54, 57, -57, 60, -60};
const int n_tilt_angles = 27;
const float SET_TILT_ANGLES[n_tilt_angles] = {0, 3, -3, 6, -6, 9, -9, 12, -12, 15, -15, 18, -18, 21, -21, 24, -24, 27, -27, 30, -30, 33, 36, 39, 42, 45, 48};
//const int n_tilt_angles = 2;
//const float SET_TILT_ANGLES[n_tilt_angles] = {60,55,50,40,45,40,35,30,25,20,15,10,5,0,-5,-10,-15,-20,-25,-30,-35,-40,-45,-50,-55,-60};
//const int n_tilt_angles = 1;
//const float SET_TILT_ANGLES[n_tilt_angles] = {55};
// Assuming a perfect counting - so no read noise, and no coincednce loss. Then we have a flat NPS and the DQE is just DQE(0)*MTF^2
// The parameters below are for a 5 gaussian fit to 300 KeV , 2.5 EPS from Ruskin et al. with DQE(0) = 0.791
const float DQE_PARAMETERS_A[1][5] = {
{-0.01516, -0.5662, -0.09731, -0.01551, 21.47}
};
const float DQE_PARAMETERS_B[1][5] = {
{0.02671, -0.02504, 0.162, 0.2831, -2.28}
};
const float DQE_PARAMETERS_C[1][5] = {
{0.01774, 0.1441, 0.1082, 0.07916, 1.372}
};
//// Rough fit for ultrascan4000
//const float DQE_PARAMETERS_A[1][5] = {
//
// {-1.356, 59.87, -0.2552, -0.00679, -0.02245 }
//
//
//};
//
//const float DQE_PARAMETERS_B[1][5] = {
//
// {-0.1878, -3.064, 0.0995, 0.2449, 0.4422}
//
//};
//
//const float DQE_PARAMETERS_C[1][5] = {
//
// {0.1207,1.457,0.264,0.006423,0.06323}
//
//};
class SimulateApp : public MyApp {
public:
bool doExpertOptions = false;
bool make_tilt_series = false;
bool need_to_allocate_projected_water = true;
bool do_beam_tilt = false;
bool use_existing_params = false;
ScatteringPotential sp;
Image* projected_water = nullptr; // waters calculated over the specified subpixel shifts and projected.
bool DoCalculation( );
void DoInteractiveUserInput( );
std::string output_filename;
float mean_defocus = 0.0f;
float defocus_1 = 0.0f;
float defocus_2 = 0.0f;
float defocus_angle = 0.0f;
float kV = 300.0f;
float spherical_aberration = 2.7f;
float objective_aperture = 100.0f;
float wavelength = 1.968697e-2; // Angstrom, default for 300kV
float wanted_pixel_size = 0.0f;
float wanted_pixel_size_sq = 0.0f;
float unscaled_pixel_size = 0.0f; // When there is a mag change
float phase_plate_thickness = 276.0; //TODO CITE Angstrom, default for pi/2 phase shift, from 2.0 g/cm^3 amorphous carbon, 300 kV and 10.7V mean inner potential as determined by electron holography (carbon atoms are ~ 8)
float astigmatism_scaling = 0.0f;
std::vector<float> image_mean;
std::vector<float> inelastic_mean;
float current_total_exposure = 0;
float pre_exposure = 0;
int size_neighborhood = 0;
int size_neighborhood_water = 0;
long make_particle_stack = 0;
int number_of_threads = 1;
bool do3d = true;
int bin3d = 0; // in addition to binning for 3d (set on input) if wanted_output_size > 0, and it is 2d, this keeps the PDB constructor from thinking it is a 3d,
float dose_per_frame = 1;
float dose_rate = 8.0; // e/physical pixel/s this should be set by the user. Changes the illumination aperature and DQE curve
float number_of_frames = 1;
double total_waters_incorporated = 0;
float average_at_cutoff[N_WATER_TERMS]; // This assumes a fixed 1/2 angstrom sampling of the hydration shell curve
float water_weight[N_WATER_TERMS];
float set_real_part_wave_function_in = 1.0f;
int minimum_padding_x_and_y = 32; // This will be changed to be the larger of this value plus the required taper region
float propagation_distance = 5.0f;
float minimum_thickness_z = 20.0f;
bool ONLY_SAVE_SUMS = true;
bool DO_NOT_RANDOMIZE_ANGLES = false;
int MODIFY_ONLY_SIGNAL = 1;
// To add error to the global alignment
float tilt_axis = 0; // degrees from Y-axis FIXME thickness calc, water padding, a few others are only valid for 0* tilt axis.
float in_plane_sigma = 2; // spread in-plane angles based on neighbors
float tilt_angle_sigma = 0.1; //;
float magnification_sigma = 0.0001; //;
float stdErr = 0.0f;
float extra_phase_shift = 0.0f;
RotationMatrix beam_tilt_xform;
float beam_tilt_x = 0.0f; //0.6f;
float beam_tilt_y = 0.0f; //-0.2f;
float beam_tilt_z_X_component;
float beam_tilt_z_Y_component;
float beam_tilt = 0.0f;
float beam_tilt_azimuth = 0.0f;
float particle_shift_x = 0.0f;
float particle_shift_y = 0.0f;
float DO_BEAM_TILT_FULL = true;
float amplitude_contrast;
float bFactor_scaling;
float min_bFactor;
FrealignParameterFile parameter_file;
cisTEMParameters input_star_file;
cisTEMParameters output_star_file;
cisTEMParameterLine parameters;
std::string parameter_file_name;
std::string input_star_file_name;
std::string output_star_file_name;
long number_preexisting_particles;
wxString preexisting_particle_file_name;
std::array<float, 17> parameter_vect;
float water_scaling = 1.0f;
float non_water_inelastic_scaling = 1.0f;
Coords coords;
StopWatch timer;
wxString wanted_symmetry = "C1";
// Intermediate images that may be useful for diagnostics.
void AddCommandLineOptions( );
bool SAVE_WATER_AND_OTHER = false;
bool SAVE_PROJECTED_WATER = false;
bool SAVE_PHASE_GRATING = false;
bool SAVE_PHASE_GRATING_DOSE_FILTERED = false;
bool SAVE_PHASE_GRATING_PLUS_WATER = false;
bool DO_CROSSHAIR = false;
bool SAVE_PROBABILITY_WAVE = false;
bool SAVE_WITH_DQE = false;
bool SAVE_WITH_NORMALIZED_DOSE = false;
bool SAVE_POISSON_PRE_NTF = false;
bool SAVE_POISSON_WITH_NTF = false;
// Add these from the command line with long option
int max_number_of_noise_particles = 0;
float noise_particle_radius_as_mutliple_of_particle_radius = 1.8;
float noise_particle_radius_randomizer_lower_bound_as_praction_of_particle_radius = -0.05;
float noise_particle_radius_randomizer_upper_bound_as_praction_of_particle_radius = 0.10;
float emulate_tilt_angle = 0.0f;
void probability_density_2d(PDB* pdb_ensemble, int time_step);
void calc_water_potential(Image* projected_water, AtomType atom_id);
void fill_water_potential(const PDB* current_specimen, Image* scattering_slab, Image* scattering_potential,
Image* inelastic_potential, Image* distance_slab, Image* water_mask_slab, Water* water_box, RotationMatrix rotate_waters,
float rotated_oZ, int* slabIDX_start, int* slabIDX_end, int iSlab);
void project(Image* image_to_project, Image* image_to_project_into, int iSlab);
void taper_edges(Image* image_to_taper, int iSlab, bool inelastic_img);
void apply_sqrt_DQE_or_NTF(Image* image_in, int iTilt_IDX, bool do_root_DQE);
//////////////////////////////////////////
////////////
// For development
//
//const float MEAN_FREE_PATH = 4000;// Angstrom, newer paper by bridgett (2018) and a couple older TODO add ref. Use to reduce total probability (https://doi.org/10.1016/j.jsb.2018.06.007)
// Note that the first two errors here are just used in matching amorphous carbon for validation. The third is used in simulations.
// The surface phase error (Wanner & Tesche 2005) quantified by holography accounts for a bias due to surface effects not included in the simple model here
float SURFACE_PHASE_ERROR = 0.497;
// The bond phase error is used to account for the remaining phase shift that is missing, due to all remaining scattering. The assumption is that amorphous water has >= the scattering due to delocalized electrons
float BOND_PHASE_ERROR = 0.0;
float BOND_SCALING_FACTOR = cistem::bond_scaling;
float MEAN_INNER_POTENTIAL = 9.09; // for 1.75 g/cm^3 amorphous carbon as in Wanner & Tesche
bool DO_PHASE_PLATE = false;
bool SHIFT_BY_CENTER_OF_MASS = true; // By default the atoms in the PDB are shifted so the center of mass is at the origin.
// CONTROL FOR DEBUGGING AND VALIDATION
bool add_mean_water_potential = false; // For comparing to published results - only applies to a 3d potential
bool DO_SINC_BLUR = false; // TODO add me back in. This is to blur the image due to intra-frame motion. Apply to projected specimen not waters
float last_shift_x = 0;
float last_shift_y = 0;
float this_shift_x = 0;
float this_shift_y = 0;
bool DO_PRINT = false;
//
bool DO_SOLVENT = true; // False to skip adding any solvent
bool CALC_WATER_NO_HOLE = false; // Makes the solvent cutoff so large that water is added on top of the protein. This is to mimick Rullgard/Vulovic
bool CALC_HOLES_ONLY = false; // This should calculate the atoms so the water is properly excluded, but not include these in the propagation. TODO checkme
bool DEBUG_POISSON = false; // Skip the poisson draw - must be true (this is gets over-ridden) if DO_PHASE_PLATE is true
bool DO_COHERENCE_ENVELOPE = true;
bool WHITEN_IMG = false; // if SAVE_REF then it will also be whitened when this option is enabled TODO checkme
bool SAVE_REF = false;
bool CORRECT_CTF = false; // only affects the image if Whiten_img is also true
bool EXPOSURE_FILTER_REF = false;
bool DO_EXPOSURE_FILTER_FINAL_IMG = false;
int DO_EXPOSURE_FILTER = 2; ///// In 2d or 3d - or 0 to turn of. 2d does not appear to produce the expected results.
bool DO_EXPOSURE_COMPLEMENT_PHASE_RANDOMIZE = false; // maintain the Energy of the protein (since no mass is actually lost) by randomizing the phases with weights that complement the exposure filter
bool DO_APPLY_DQE = true;
int CAMERA_MODEL = 0;
int wanted_output_size = -1;
float wgt = 0.0f;
float bf = 0.0f;
bool water_shell_only = true;
bool is_alpha_fold_prediction = false;
///////////
/////////////////////////////////////////
private:
};
// Optional command-line stuff
void SimulateApp::AddCommandLineOptions( ) {
// TODO consider short vs long switches.
// Options for saving diagnostic images
command_line_parser.AddLongSwitch("save-water-and-other", "Save image of water and scattering.");
command_line_parser.AddLongSwitch("save-projected-water", "Save image projected water atoms with sub-pixel offsets.");
command_line_parser.AddLongSwitch("save-phase-grating", "Save phase grating of non-solvent atoms.");
command_line_parser.AddLongSwitch("do-cross-hair", "Trim to output size and abort after phase calculation.");
command_line_parser.AddLongSwitch("save-phase-grating-dose-filtered", "Save phase grating of non-solvent atoms + exposure filter.");
command_line_parser.AddLongSwitch("save-phase-grating-plus-water", "Save phase grating of non-solvent atoms plus water atoms.");
command_line_parser.AddLongSwitch("save-with-dqe", "Save |wavefunction|^2 + DQE effects.");
// These shouldn't be needed anymore, no NTF applied.
command_line_parser.AddLongSwitch("save-poisson-pre-ntf", "Save image of water and scattering.");
command_line_parser.AddLongSwitch("save-poisson-post-ntf", "Save image of water and scattering.");
command_line_parser.AddLongSwitch("save-detector-wavefunction", "Save the detector wave function directly? Skip Poisson draw, default is no"); // Skip the poisson draw - must be true (this is gets over-ridden) if DO_PHASE_PLATE is true
// Option to skip centering by mass
command_line_parser.AddLongSwitch("skip-centering-by-mass", "Skip centering by mass");
// Options for computation (TODO add gpu flag here)
command_line_parser.AddOption("j", "", "Desired number of threads. Overrides interactive user input. Is overriden by env var OMP_NUM_THREADS", wxCMD_LINE_VAL_NUMBER);
// Options related to calibration
command_line_parser.AddOption("", "surface-phase-error", "Surface phase error (radian) From Wanner & Tesche, default is 0.497", wxCMD_LINE_VAL_DOUBLE);
command_line_parser.AddOption("", "bond-phase-error", "Bond phase error (radian) Calibrated from experiment, may need to be updated, default is 0.09", wxCMD_LINE_VAL_DOUBLE);
// Note that the first two errors here are just used in matching amorphous carbon for validation. The third is used in simulations.
// The surface phase error (Wanner & Tesche 2005) quantified by holography accounts for a bias due to surface effects not included in the simple model here
// The bond phase error is used to account for the remaining phase shift that is missing, due to all remaining scattering. The assumption is that amorphous water has >= the scattering due to delocalized electrons
// To account for the bond phase error in practice, a small scaling factor is applied to the atomic potentials
command_line_parser.AddOption("", "bond-scaling", "Compensate for bond phase error by scaling the interaction potential. Calibrated from experiment, may need to be updated. default is 1.065", wxCMD_LINE_VAL_DOUBLE);
command_line_parser.AddOption("", "carbon-mip", "Amorphous carbon MIP From Wanner & Tesche for 1.75 g/cm^3 carbon. default is 9.09", wxCMD_LINE_VAL_DOUBLE); // for 1.75 g/cm^3 amorphous carbon as in Wanner & Tesche
command_line_parser.AddLongSwitch("do-phase-plate", "do a phase plate simulation? default is no");
command_line_parser.AddOption("", "non-water-inelastic-scaling", "Scale the apparent atomic number for water from its default of 7.35", wxCMD_LINE_VAL_DOUBLE);
command_line_parser.AddOption("", "wanted-symmetry", "Desired symmetry.", wxCMD_LINE_VAL_STRING);
// CONTROL FOR DEBUGGING AND VALIDATION
command_line_parser.AddLongSwitch("add-constant-background", "Add a constant potential for the mean water? default is no"); // For comparing to published results - only applies to a 3d potential
command_line_parser.AddLongSwitch("do-sinc-blur", "Intra-frame sinc blur on protein? Currently this does nothing, default is no"); // TODO add me back in. This is to blur the image due to intra-frame motion. Apply to projected specimen not waters
command_line_parser.AddLongSwitch("print-extra-info", "Print out a bunch of extra diagnostic information? default is no");
command_line_parser.AddLongSwitch("skip-solvent", "Skip adding potential, constant or atom-wise for solvent? default is no"); // False to skip adding any solvent
command_line_parser.AddLongSwitch("calc-water-no-hole", "Ignore protein envelope and add solvent everywhere? default is no"); // Makes the solvent cutoff so large that water is added on top of the protein. This is to mimick Rullgard/Vulovic
command_line_parser.AddLongSwitch("calc-holes-only", "Use protein envelope, but don't add it in (holes only)? default is no"); // This should calculate the atoms so the water is properly excluded, but not include these in the propagation. TODO checkme
// command_line_parser.AddLongSwitch("save-detector-wavefunction", "Save the detector wave function directly? Skip Poisson draw, default is no"); // Skip the poisson draw - must be true (this is gets over-ridden) if DO_PHASE_PLATE is true
command_line_parser.AddLongSwitch("skip-random-angles", "Skip randomizing angles in a particle stack? default is no"); //
// command_line_parser.AddLongSwitch("only-modify-signal-3d", "When applying the cummulative exposure filter to the 3d reference, only modify the signal (reduce exposure filter by 1 - (1-x)/(1+x)");
command_line_parser.AddOption("", "only-modify-signal-3d", "When applying the cummulative exposure filter to the 3d reference, only modify the signal (reduce exposure filter by 1 - (1-x)/(1+x)", wxCMD_LINE_VAL_NUMBER);
command_line_parser.AddLongSwitch("skip-coherence-envelope", "Apply spatial coherence envelope? default is no"); // These do nothing~! FIXME
command_line_parser.AddOption("", "max_number_of_noise_particles", "Maximum number of neighboring noise particles when simulating an image stack. Default is 0", wxCMD_LINE_VAL_NUMBER);
command_line_parser.AddOption("", "noise_particle_radius_as_mutliple_of_particle_radius", "default is 1.8", wxCMD_LINE_VAL_DOUBLE);
command_line_parser.AddOption("", "noise_particle_radius_randomizer_lower_bound_as_praction_of_particle_radius", "default is -0.05", wxCMD_LINE_VAL_DOUBLE);
command_line_parser.AddOption("", "noise_particle_radius_randomizer_upper_bound_as_praction_of_particle_radius", "default is 0.10", wxCMD_LINE_VAL_DOUBLE);
command_line_parser.AddOption("", "emulate_tilt_angle", "default is 0.0 degrees around Y axis", wxCMD_LINE_VAL_DOUBLE);
command_line_parser.AddLongSwitch("whiten-output", "Whiten the image? default is no"); // if SAVE_REF then it will also be whitened when this option is enabled TODO checkme
command_line_parser.AddLongSwitch("do-perfect-reference", "Save a perfect reference image with no noise in addition to the image. default is no");
command_line_parser.AddLongSwitch("apply-ctf", "Multiply by average CTF? default is no"); // only affects the image if Whiten_img is also true
command_line_parser.AddLongSwitch("apply-exposure-filter-ref", "Exposure filter the perfect reference? default is no");
command_line_parser.AddLongSwitch("apply-exposure-filter-img", "Exposure filter the final movie/image? default is no");
command_line_parser.AddLongSwitch("skip-radiation-damage", "Do NOT apply the 2d exposure filter to the scattering potential prior to adding the solvent. default is no"); ///// In 2d or 3d - or 0 to turn of. 2d does not appear to produce the expected results.
command_line_parser.AddLongSwitch("maintain-power", "Maintain power with random phases in exposure filter? default is no"); // maintain the Energy of the protein (since no mass is actually lost) by randomizing the phases with weights that complement the exposure filter
command_line_parser.AddLongSwitch("skip-dqe", "Do NOT apply the DQE? depends on camera model default is no and applies K2 fixme");
command_line_parser.AddLongSwitch("skip-tilted-propagation", "Apply the phase shift due to beam tilt, but don't actually do the inclined wave propagation. default no");
command_line_parser.AddLongSwitch("save-frames", "The default is to save the integrated frames. This option overrides this behavior.");
command_line_parser.AddOption("", "wgt", "Maximum number of neighboring noise particles when simulating an image stack. Default is 0", wxCMD_LINE_VAL_DOUBLE);
command_line_parser.AddOption("", "bf", "Maximum number of neighboring noise particles when simulating an image stack. Default is 0", wxCMD_LINE_VAL_DOUBLE);
command_line_parser.AddLongSwitch("disable-water-shell-only", "normally when adding constant background to a 3d, taper off 4 Ang into the water");
command_line_parser.AddLongSwitch("is-alpha-fold-prediction", "Is this a alpha-fold prediction? If so, convert the confidence score stored in the bfactor column to a bfactor. default is no");
// command_line_parser.AddOption("j","","Desired number of threads. Overrides interactive user input. Is overriden by env var OMP_NUM_THREADS",wxCMD_LINE_VAL_NUMBER);
}
IMPLEMENT_APP(SimulateApp);
// override the DoInteractiveUserInput
void SimulateApp::DoInteractiveUserInput( ) {
bool add_more_pdbs = true;
bool supply_origin = false;
int iPDB = 0;
int iOrigin;
int iParticleCopy;
double temp_double;
long temp_long;
//////////
/////////// Check the command line options, could this be done at the top?
SAVE_WATER_AND_OTHER = command_line_parser.FoundSwitch("save-water-and-other");
SAVE_PROJECTED_WATER = command_line_parser.FoundSwitch("save-projected-water");
SAVE_PHASE_GRATING = command_line_parser.FoundSwitch("save-phase-grating");
DO_CROSSHAIR = command_line_parser.FoundSwitch("do-cross-hair"); // implies skip-solvent
SAVE_PHASE_GRATING_DOSE_FILTERED = command_line_parser.FoundSwitch("save-phase-grating-dose-filtered");
SAVE_PHASE_GRATING_PLUS_WATER = command_line_parser.FoundSwitch("save-phase-grating-plus-water");
SAVE_WITH_DQE = command_line_parser.FoundSwitch("save-with-dqe");
SAVE_POISSON_PRE_NTF = command_line_parser.FoundSwitch("save-poisson-pre-ntf");
SAVE_POISSON_WITH_NTF = command_line_parser.FoundSwitch("save-poisson-post-ntf");
DEBUG_POISSON = command_line_parser.Found("save-detector-wavefunction");
if ( command_line_parser.Found("j", &temp_long) ) {
number_of_threads = (int)temp_long;
}
if ( command_line_parser.Found("surface-phase-error", &temp_double) ) {
SURFACE_PHASE_ERROR = (float)temp_double;
}
if ( command_line_parser.Found("bond-phase-error", &temp_double) ) {
BOND_PHASE_ERROR = (float)temp_double;
}
// Note that the first two errors here are just used in matching amorphous carbon for validation. The third is used in simulations.
// The surface phase error (Wanner & Tesche 2005) quantified by holography accounts for a bias due to surface effects not included in the simple model here
// The bond phase error is used to account for the remaining phase shift that is missing, due to all remaining scattering. The assumption is that amorphous water has >= the scattering due to delocalized electrons
// To account for the bond phase error in practice, a small scaling factor is applied to the atomic potentials
if ( command_line_parser.Found("bond-scaling", &temp_double) ) {
BOND_SCALING_FACTOR = (float)temp_double;
}
if ( command_line_parser.Found("carbon-mip", &temp_double) ) {
MEAN_INNER_POTENTIAL = (float)temp_double;
}
if ( command_line_parser.FoundSwitch("do-phase-plate") )
DO_PHASE_PLATE = true;
if ( command_line_parser.Found("non-water-inelastic-scaling", &temp_double) ) {
non_water_inelastic_scaling = (float)temp_double;
}
if ( command_line_parser.FoundSwitch("skip-centering-by-mass") )
SHIFT_BY_CENTER_OF_MASS = false;
// CONTROL FOR DEBUGGING AND VALIDATION
add_mean_water_potential = command_line_parser.Found("add-constant-background");
DO_SINC_BLUR = command_line_parser.Found("do-sinc-blur");
DO_PRINT = command_line_parser.Found("print-extra-info");
if ( command_line_parser.Found("skip-solvent") )
DO_SOLVENT = false;
CALC_WATER_NO_HOLE = command_line_parser.Found("calc-water-no-hole");
CALC_HOLES_ONLY = command_line_parser.Found("calc-holes-only");
if ( command_line_parser.Found("skip-coherence-envelope") )
DO_COHERENCE_ENVELOPE = false;
WHITEN_IMG = command_line_parser.Found("whiten-output");
SAVE_REF = command_line_parser.Found("do-perfect-reference");
CORRECT_CTF = command_line_parser.Found("apply-ctf");
EXPOSURE_FILTER_REF = command_line_parser.Found("apply-exposure-filter-ref");
DO_EXPOSURE_FILTER_FINAL_IMG = command_line_parser.Found("apply-exposure-filter-img");
if ( command_line_parser.Found("skip-radiation-damage") )
DO_EXPOSURE_FILTER = false;
DO_EXPOSURE_COMPLEMENT_PHASE_RANDOMIZE = command_line_parser.Found("maintain-power");
if ( command_line_parser.Found("skip-dqe") )
DO_APPLY_DQE = false;
if ( command_line_parser.Found("skip-tilted-propagation") )
DO_BEAM_TILT_FULL = false;
if ( command_line_parser.Found("save-frames") )
ONLY_SAVE_SUMS = false;
if ( command_line_parser.Found("skip-random-angles") )
DO_NOT_RANDOMIZE_ANGLES = true;
if ( command_line_parser.Found("only-modify-signal-3d", &temp_long) ) {
MODIFY_ONLY_SIGNAL = (int)temp_long;
}
if ( command_line_parser.Found("max_number_of_noise_particles", &temp_long) ) {
max_number_of_noise_particles = (int)temp_long;
}
if ( command_line_parser.Found("noise_particle_radius_as_mutliple_of_particle_radius", &temp_double) ) {
noise_particle_radius_as_mutliple_of_particle_radius = (float)temp_double;
}
if ( command_line_parser.Found("noise_particle_radius_randomizer_lower_bound_as_praction_of_particle_radius", &temp_double) ) {
noise_particle_radius_randomizer_lower_bound_as_praction_of_particle_radius = (float)temp_double;
}
if ( command_line_parser.Found("noise_particle_radius_randomizer_upper_bound_as_praction_of_particle_radius", &temp_double) ) {
noise_particle_radius_randomizer_upper_bound_as_praction_of_particle_radius = (float)temp_double;
}
if ( command_line_parser.Found("emulate_tilt_angle", &temp_double) ) {
emulate_tilt_angle = (float)temp_double;
}
wxString temp_string;
if ( command_line_parser.Found("wanted-symmetry", &temp_string) ) {
wanted_symmetry = temp_string;
}
if ( command_line_parser.Found("wgt", &temp_double) ) {
wgt = (float)temp_double;
}
if ( command_line_parser.Found("bf", &temp_double) ) {
bf = (float)temp_double;
}
if ( command_line_parser.Found("disable-water-shell-only") )
water_shell_only = false;
if ( command_line_parser.Found("is-alpha-fold-prediction") )
is_alpha_fold_prediction = true;
if ( DO_CROSSHAIR )
DO_SOLVENT = false;
UserInput* my_input = new UserInput("Simulator", 0.25);
output_filename = my_input->GetFilenameFromUser("Output filename", "just the base, no extension, will be mrc", "test_tilt.mrc", false);
// FIXME the range is way too big for 3d but needed for the fixed size hack.
do3d = my_input->GetYesNoFromUser("Make a 3d scattering potential?", "just potential if > 0, input is the wanted cubic size", "yes");
wanted_output_size = my_input->GetIntFromUser("Wanted output size", "Size along all dimensions", "16", 1, 8096);
// Limit the total 3d size allowed. This should probably depend on the pixel size, because of the oversampling
if ( do3d && wanted_output_size > MAX_3D_SIZE ) {
wxPrintf("WARNING: wanted 3d size is too big, max size is %d\n", MAX_3D_SIZE);
exit(-1);
}
number_of_threads = my_input->GetIntFromUser("Number of threads", "Max is number of tilts", "1", 1);
bool require_pdb_to_exist = true;
int runaway_counter = 0;
while ( add_more_pdbs && runaway_counter < 10000 ) {
sp.pdb_file_names.push_back(my_input->GetFilenameFromUser("PDB file name", "an input PDB", "my_atom_coords.pdb", require_pdb_to_exist));
add_more_pdbs = my_input->GetYesNoFromUser("Add another type of particle?", "Add another pdb to create additional features in the ensemble", "no");
runaway_counter++;
if ( runaway_counter > 1000 ) {
wxPrintf("ERROR: runaway_counter > 1000");
exit(-1);
}
}
// Allow bond scaling factor (constants/cistem::bond_scaling) to be overrode for experiment.
sp.SetImagingParameters(wanted_pixel_size, kV, BOND_SCALING_FACTOR);
wanted_pixel_size = my_input->GetFloatFromUser("Output pixel size (Angstroms)", "Output size for the final projections", "1.0", 0.01, MAX_PIXEL_SIZE);
bFactor_scaling = my_input->GetFloatFromUser("Linear scaling of per atom bFactor", "0 off, 1 use as is", "0", 0, 10000);
min_bFactor = my_input->GetFloatFromUser("Per atom (xtal) bFactor added to all atoms", "accounts for all quote[unavoidable] experimental error", "15.0", 0.0f, 10000);
if ( do3d ) {
bool found_the_best_binning = false;
// Check to make sure the sampling is sufficient, if not, oversample and bin at the end.
if ( wanted_pixel_size > 1.5 ) {
if ( wanted_output_size * 4 < MAX_3D_SIZE ) {
wxPrintf("\nOversampling your 3d by a factor of 4 for calculation.\n");
wanted_pixel_size /= 4.f;
bin3d = 4;
found_the_best_binning = true;
}
}
if ( wanted_pixel_size > 0.75 && ! found_the_best_binning ) {
if ( wanted_output_size * 2 < MAX_3D_SIZE ) {
wxPrintf("\nOversampling your 3d by a factor of 2 for calculation.\n");
wanted_pixel_size /= 2.f;
bin3d = 2;
found_the_best_binning = true;
}
}
if ( ! found_the_best_binning ) {
bin3d = 1;
}
dose_per_frame = my_input->GetFloatFromUser("electrons/Ang^2 in a frame at the specimen", "", "1.0", 0.05, 20.0);
number_of_frames = my_input->GetFloatFromUser("number of frames per movie (micrograph or tilt)", "", "30", 1.0, 1000.0);
}
else {
// First we check on whether there is an input set of parameters, as this changes the most options following
use_existing_params = my_input->GetYesNoFromUser("Use an existing set of orientations", "yes no", "no");
if ( use_existing_params ) {
preexisting_particle_file_name = my_input->GetFilenameFromUser("cisTEM star file name", "an input star file to match reconstruction", "myparams.star", true);
int default_number_parameters = 1;
if ( DoesFileExist(preexisting_particle_file_name) ) {
input_star_file.ReadFromcisTEMStarFile(preexisting_particle_file_name);
number_preexisting_particles = input_star_file.ReturnNumberofLines( );
default_number_parameters = number_preexisting_particles;
wxPrintf("\nFound %ld particles in the input star file\n", number_preexisting_particles);
}
else {
SendErrorAndCrash(wxString::Format("Error: Input star file %s not found\n", preexisting_particle_file_name));
}
// The following doesn't work if there is a .dff file
// std::string wanted_default = std::to_string(default_number_parameters);
// number_preexisting_particles = my_input->GetIntFromUser("Limit the number of particles used from the star file.", "0: use all", wanted_default.c_str());
}
// image type :
// The default option is to simulate a full micrograph. Ask if instead we want a particle stack or a tilt-series
make_particle_stack = (long)my_input->GetFloatFromUser("Create a particle stack rather than a micrograph?", "Number of particles at random orientations, 0 for just an image", "1", 0, 1e7);
make_tilt_series = my_input->GetYesNoFromUser("Create a tilt-series as well?", "Note that this is mutually exclusive with creating a particle stack", "no");
if ( make_tilt_series && make_particle_stack ) {
SendErrorAndCrash("Error: You can't create a tilt-series and a particle stack at the same time");
}
// This is currently undefined behavior
if ( make_tilt_series && use_existing_params ) {
SendErrorAndCrash("Error: You can't create a tilt-series and use an existing set of orientations at the same time");
}
bool wanted_parameters_are_present = false;
if ( use_existing_params ) {
// If we are using existing parameters, we enforce the following are specified.
if ( input_star_file.parameters_that_were_read.defocus_1 &&
input_star_file.parameters_that_were_read.defocus_2 &&
input_star_file.parameters_that_were_read.defocus_angle &&
input_star_file.parameters_that_were_read.phase_shift &&
input_star_file.parameters_that_were_read.pixel_size &&
input_star_file.parameters_that_were_read.pre_exposure &&
input_star_file.parameters_that_were_read.total_exposure &&
input_star_file.parameters_that_were_read.microscope_spherical_aberration_mm &&
input_star_file.parameters_that_were_read.microscope_voltage_kv &&
input_star_file.parameters_that_were_read.beam_tilt_x &&
input_star_file.parameters_that_were_read.beam_tilt_y ) {
wanted_parameters_are_present = true;
kV = input_star_file.ReturnMicroscopekV(0);
spherical_aberration = input_star_file.ReturnMicroscopeCs(0);
float max_exposure, mean_phase_shift, mean_beam_tilt_x(0.f), mean_beam_tilt_y(0.f);
mean_defocus = 0.f;
mean_phase_shift = 0.f;
number_of_frames = 1;
// Assuming this is constant for all particles and frames FIXME
dose_per_frame = input_star_file.ReturnTotalExposure(0) - input_star_file.ReturnPreExposure(0);
// tmp override for testing FIXME
dose_per_frame = 1;
current_total_exposure = 1;
pre_exposure = 1;
for ( int counter = 0; counter < number_preexisting_particles; counter++ ) {
mean_defocus += 0.5f * (input_star_file.ReturnDefocus1(counter) + input_star_file.ReturnDefocus2(counter));
number_of_frames = std::max(number_of_frames, float(input_star_file.ReturnParticleGroup(counter))); // FIXME, why is number of frames a float, prob a bad idea
mean_beam_tilt_x += input_star_file.ReturnBeamTiltX(counter);
mean_beam_tilt_y += input_star_file.ReturnBeamTiltY(counter);
}
mean_defocus /= number_preexisting_particles;
beam_tilt_x = mean_beam_tilt_x / number_preexisting_particles;
beam_tilt_y = mean_beam_tilt_y / number_preexisting_particles;
dose_rate = 3.0f; // FIXME this is not currently in the parameter file
output_star_file = input_star_file; // we may change some of the parameters or save frames
}
else {
SendInfo("Warning: You must specify all parameters in the input star file\nIgnoring the input star file.\nQuit, or continue by entering manual parameters.\n");
}
}
if ( ! wanted_parameters_are_present ) {
// Either because the user did not ask for them (no-preexisting params) or because they were not present in the input star file
mean_defocus = my_input->GetFloatFromUser("wanted mean_defocus (Angstroms)", "Out", "700", 0, 120000);
extra_phase_shift = my_input->GetFloatFromUser("wanted additional phase shift x * PI radians, i.e. 0.5 for PI/2 shift.", "", "0.0", -2.0, 2.0);
dose_per_frame = my_input->GetFloatFromUser("electrons/Ang^2 in a frame at the specimen", "", "1.0", 0.05, 20.0);
dose_rate = my_input->GetFloatFromUser("electrons/Pixel/sec", "Affects coherence but not coincidence loss", "3.0", 0.001, 200.0);
number_of_frames = my_input->GetFloatFromUser("number of frames per movie (micrograph or tilt)", "", "30", 1.0, 1000.0);
}
}
doExpertOptions = my_input->GetYesNoFromUser("Set expert options?", "", "no");
wanted_pixel_size_sq = powf(wanted_pixel_size, 2);
if ( doExpertOptions ) {
water_scaling = my_input->GetFloatFromUser("Linear scaling water intensity", "0 off, 1 use as is", "1", 0, 1);
astigmatism_scaling = my_input->GetFloatFromUser("fraction of the mean_defocus to add as astigmatism", "0 off, 1 use as is", "0.0", 0, 0.5);
// FIXME with star file reading these are just silently ignored, but should probably n ot even be asked for
if ( use_existing_params ) {
my_input->GetFloatFromUser("IGNORED with input star: Accelrating volatage (kV)", "Default is 300", "300.0", 80.0, 1000.0f); // Calculations are not valid outside this range
my_input->GetFloatFromUser("IGNORED iwth input star: Spherical aberration constant in millimeters", "", "2.7", 0.0, 5.0);
}
else {
kV = my_input->GetFloatFromUser("Accelrating volatage (kV)", "Default is 300", "300.0", 80.0, 1000.0f); // Calculations are not valid outside this range
spherical_aberration = my_input->GetFloatFromUser("Spherical aberration constant in millimeters", "", "2.7", 0.0, 5.0);
}
objective_aperture = my_input->GetFloatFromUser("Objective aperture diameter in micron", "", "100.0", 0.0, 1000.0);
stdErr = my_input->GetFloatFromUser("Std deviation of error to use in shifts, astigmatism, rotations etc.", "", "0.0", 0.0, 100.0);
pre_exposure = my_input->GetFloatFromUser("Pre-exposure in electron/A^2", "use for testing exposure filter", "0", 0.0);
// This minimum padding should probably depend on the mean_defocus. FIXME
minimum_padding_x_and_y = my_input->GetIntFromUser("minimum padding of images with solvent", "", "32", 0, 4096);
minimum_thickness_z = my_input->GetIntFromUser("minimum thickness in Z", "", "10", 2, 10000);
propagation_distance = my_input->GetFloatFromUser("Propagation distance in angstrom", "Also used as minimum thickness", "5", -1e6, 1e6);
if ( fabsf(propagation_distance) > minimum_thickness_z ) {
minimum_thickness_z = fabsf(propagation_distance);
if ( DO_PRINT ) {
wxPrintf("min thickness was less than propagation distance, so setting it there\n");
}
}
// To add error to the global alignment
tilt_axis = my_input->GetFloatFromUser("Rotation of tilt-axis from Y (degrees)", "", "0.0", 0, 180); // TODO does this apply everywhere it should
in_plane_sigma = my_input->GetFloatFromUser("Standard deviation on angles in plane (degrees)", "", "2", 0, 100); // spread in-plane angles based on neighbors
tilt_angle_sigma = my_input->GetFloatFromUser("Standard deviation on tilt-angles (degrees)", "", "0.1", 0, 10);
; //;
magnification_sigma = my_input->GetFloatFromUser("Standard deviation on magnification (fraction)", "", "0.0001", 0, 1); //;
if ( ! use_existing_params ) {
beam_tilt_x = my_input->GetFloatFromUser("Beam-tilt in X (milli radian)", "", "0.0", -300, 300); //0.6f;
beam_tilt_y = my_input->GetFloatFromUser("Beam-tilt in Y (milli radian)", "", "0.0", -300, 300); //-0.2f;
particle_shift_x = my_input->GetFloatFromUser("Beam-tilt particle shift in X (Angstrom)", "", "0.0", -100, 100);
particle_shift_y = my_input->GetFloatFromUser("Beam-tilt particle shift in Y (Angstrom)", "", "0.0", -100, 100);
}
set_real_part_wave_function_in = sqrt(dose_per_frame * wanted_pixel_size_sq);
}
else {
water_scaling = 1.0f;
astigmatism_scaling = 0.0f;
objective_aperture = 100.0f;
spherical_aberration = 2.7f;
stdErr = 0.0f;
}
wavelength = ReturnWavelenthInAngstroms(kV);
output_star_file_name = output_filename + ".star";
if ( DO_PHASE_PLATE ) {
if ( SURFACE_PHASE_ERROR < 0 ) {
if ( DO_PRINT ) {
wxPrintf("SURFACE_PHASE_ERROR < 0, subbing min thickness for phase plate thickness\n");
}
phase_plate_thickness = minimum_thickness_z;
}
else {
phase_plate_thickness = (pi_v<float> / 2.0f + SURFACE_PHASE_ERROR + BOND_PHASE_ERROR) / (MEAN_INNER_POTENTIAL / (kV * 1000) * (511 + kV) / (2 * 511 + kV) * (2 * pi_v<float> / (wavelength * 1e-10))) * 1e10;
}
if ( DO_PRINT ) {
wxPrintf("With a mean inner potential of %2.2fV a thickness of %2.2f ang is needed for a pi/2 phase shift \n", MEAN_INNER_POTENTIAL, phase_plate_thickness);
}
if ( DO_PRINT ) {
wxPrintf("Phase shift params %f %f %f\n", BOND_PHASE_ERROR, BOND_SCALING_FACTOR, SURFACE_PHASE_ERROR);
}
}
// The third term is a rough estiamte to ensure any delocalized info from particles is retained. It should probably also consider the stdErr
minimum_padding_x_and_y = myroundint(N_TAPERS * TAPERWIDTH * wanted_pixel_size + (float)minimum_padding_x_and_y + mean_defocus / 100.0f);
non_water_inelastic_scaling *= inelastic_scalar_water;
// Various overrides based on input flags and command line options.
if ( do3d ) {
// We only want one slab for a 3d simulation
minimum_thickness_z = wanted_output_size * wanted_pixel_size;
propagation_distance = -minimum_thickness_z;
// We need to be sure there is SOME exposure. Note for 3d sims, the exposure applied is an integrated
// exposure.
if ( number_of_frames - dose_per_frame < 1 ) {
number_of_frames++;
}
}
delete my_input;
}
// overide the do calculation method which will be what is actually run..
bool SimulateApp::DoCalculation( ) {
// TODO is this the best place to put this?
current_total_exposure = pre_exposure;
if ( CORRECT_CTF && use_existing_params ) {
wxPrintf("I did not set up ctf correction and the use of existing parameters. FIXME\n");
exit(-1);
}
sp.SetImagingParameters(wanted_pixel_size, kV);
if ( make_particle_stack > 0 ) {
sp.InitPdbEnsemble(SHIFT_BY_CENTER_OF_MASS, minimum_padding_x_and_y, minimum_thickness_z,
max_number_of_noise_particles,
noise_particle_radius_as_mutliple_of_particle_radius,
noise_particle_radius_randomizer_lower_bound_as_praction_of_particle_radius,
noise_particle_radius_randomizer_upper_bound_as_praction_of_particle_radius,
emulate_tilt_angle, is_alpha_fold_prediction, input_star_file, use_existing_params);
}
else {
// Over-ride the max number of noise particles
max_number_of_noise_particles = 0;
sp.InitPdbEnsemble(SHIFT_BY_CENTER_OF_MASS, minimum_padding_x_and_y, minimum_thickness_z,
max_number_of_noise_particles,
noise_particle_radius_as_mutliple_of_particle_radius,
noise_particle_radius_randomizer_lower_bound_as_praction_of_particle_radius,
noise_particle_radius_randomizer_upper_bound_as_praction_of_particle_radius,
emulate_tilt_angle, is_alpha_fold_prediction, input_star_file, use_existing_params);
}
if ( DO_PRINT ) {
wxPrintf("\nThere are %ld non-water atoms in the specimen.\n", sp.ReturnTotalNumberOfNonWaterAtoms( ));
}
if ( DO_PRINT ) {
wxPrintf("\nCurrent number of PDBs %d\n", sp.pdb_file_names.size( ));
}
// FIXME add time steps.
int time_step = 0;
this->probability_density_2d(sp.pdb_ensemble.data( ), time_step);
if ( DO_PRINT ) {
wxPrintf("\nFinished pre seg fault\n");
}
// It gives a segfault at the end either way.
// pdb_ensemble[0].Close();
return true;
}
/*
I've moved the wanted originz and euler angles as REMARK 351 in the PDB so that environments may be "easily" created in chimera.
It makes more sense then to intialize the trajectories in the call to PDB::init
Leave this in until convinced it works ok.
*/
void SimulateApp::probability_density_2d(PDB* pdb_ensemble, int time_step) {
RandomNumberGenerator my_rand(pi_v<float>);
bool SCALE_DEFOCUS_TO_MATCH_300 = true;
float scale_defocus = 1.0f;
// TODO Set even range in z to avoid large zero areas
// TODO Set a check on the solvent fraction and scaling and report if it is unreasonable. Define reasonable
// TODO Set a check on the range of values, report if mean_defocus tolerance is too small (should all be positive)
// Image img;
// img.QuickAndDirtyReadSlice("/groups/grigorieff/home/himesb/cisTEM_2/cisTEM/trunk/gpu/include/oval_full.mrc",1);
// img.QuickAndDirtyWriteSlice("/groups/grigorieff/home/himesb/tmp/noshift.mrc",1,1,false);
// img.PhaseShift(1.5,3.5,0.0);
// img.QuickAndDirtyWriteSlice("/groups/grigorieff/home/himesb/tmp/withShift.mrc",1,1,false);
// exit(-1);
long current_atom;
long nOutOfBounds = 0;
long iTilt_IDX;
int iSlab = 0;
int current_3D_slice_to_save = 0;
std::vector<float> shift_x;
std::vector<float> shift_y;
std::vector<float> shift_z;
std::vector<float> mag_diff;
float euler1(0), euler2(0), euler3(0);
// CTF parameters: There should be an associated variablility with tilt angle TODO and get rid of extra parameters
float wanted_acceleration_voltage = this->kV; // keV
float wanted_spherical_aberration = this->spherical_aberration; // mm
float wanted_defocus_1_in_angstroms = this->mean_defocus; // A
float wanted_defocus_2_in_angstroms = this->mean_defocus; //A
float wanted_astigmatism_azimuth = 0.0; // degrees
float astigmatism_angle_randomizer = 0.0; //
float defocus_randomizer = 0.0;
float wanted_additional_phase_shift_in_radians = this->extra_phase_shift * PI;
std::vector<float> propagator_distance; // in Angstom, <= mean_defocus tolerance.
float defocus_offset = 0;
if ( DO_PRINT ) {
wxPrintf("Using extra phase shift of %f radians\n", wanted_additional_phase_shift_in_radians);
}
int frame_lines;
if ( ONLY_SAVE_SUMS ) {
frame_lines = 1;
}
else {
frame_lines = number_of_frames;
}
if ( use_existing_params && make_particle_stack == 0 ) {
// We are making a micrograph with multiple particles, so scale by that number.
frame_lines *= number_preexisting_particles;
}
parameter_vect[0] = 1; // idx
parameter_vect[1] = 0; // psi
parameter_vect[2] = 0; // theta
parameter_vect[3] = 0; // phi
parameter_vect[4] = 0; // shx
parameter_vect[5] = 0; // shy
parameter_vect[6] = 0; // mag
parameter_vect[7] = 1; // include
parameter_vect[8] = wanted_defocus_1_in_angstroms; //
parameter_vect[9] = wanted_defocus_2_in_angstroms; //
parameter_vect[10] = wanted_astigmatism_azimuth; //
parameter_vect[11] = wanted_additional_phase_shift_in_radians;
parameter_vect[12] = 100; // Occupancy
parameter_vect[13] = -1000; // LogP
parameter_vect[14] = 10; //Sigma
parameter_vect[15] = 10; //Score
parameter_vect[16] = 0; // Change
// Keep a copy of the unscaled pixel size to handle magnification changes.
this->unscaled_pixel_size = this->wanted_pixel_size;
// Keep a copy of the unscaled pixel size to handle magnification changes.
this->unscaled_pixel_size = this->wanted_pixel_size;
// Setup the output star file
if ( this->make_tilt_series ) {
output_star_file.PreallocateMemoryAndBlank(n_tilt_angles * frame_lines);
}
else if ( this->make_particle_stack > 0 ) {
output_star_file.PreallocateMemoryAndBlank(make_particle_stack * frame_lines);
}
else {
// if we are using pre-existing parameters and not a particle stack, frame_lines is already scaled by the number of particles.
output_star_file.PreallocateMemoryAndBlank(frame_lines);
}
output_star_file.parameters_to_write.SetAllToTrue( );
output_star_file.parameters_to_write.image_is_active = false;
output_star_file.parameters_to_write.original_image_filename = false;
output_star_file.parameters_to_write.reference_3d_filename = false;
output_star_file.parameters_to_write.stack_filename = false;
// TODO add a method to set defaults and manage this paramter line
parameters.position_in_stack = 0;
parameters.image_is_active = 0;
parameters.psi = 0.0f;
parameters.theta = 0.0f;
parameters.phi = 0.0f;
parameters.x_shift = 0.0f;
parameters.y_shift = 0.0f;
parameters.defocus_1 = wanted_defocus_1_in_angstroms;
parameters.defocus_2 = wanted_defocus_2_in_angstroms;
parameters.defocus_angle = wanted_astigmatism_azimuth;
parameters.phase_shift = wanted_additional_phase_shift_in_radians;
parameters.occupancy = 100.0f;
parameters.logp = -1000.0f;
parameters.sigma = 10.0f;
parameters.score = 10.0f;
parameters.score_change = 0.0f;
parameters.pixel_size = wanted_pixel_size;
parameters.microscope_voltage_kv = wanted_acceleration_voltage;
parameters.microscope_spherical_aberration_mm = wanted_spherical_aberration;
parameters.amplitude_contrast = 0.0f;
parameters.beam_tilt_x = beam_tilt_x;
parameters.beam_tilt_y = beam_tilt_y;
parameters.image_shift_x = particle_shift_x;
parameters.image_shift_y = particle_shift_y;
parameters.stack_filename = output_filename;
parameters.original_image_filename = wxEmptyString;
parameters.reference_3d_filename = wxEmptyString;
parameters.best_2d_class = 0;
parameters.beam_tilt_group = 0;
parameters.particle_group = 0;
parameters.pre_exposure = 0.0f;
parameters.total_exposure = 0.0f;
// Do this after intializing the parameters which should be stored in millirad
if ( beam_tilt_x != 0.0 || beam_tilt_y != 0.0 ) {
do_beam_tilt = true;
beam_tilt_azimuth = atan2f(beam_tilt_y, beam_tilt_x);
beam_tilt = sqrtf(powf(beam_tilt_x, 2) + powf(beam_tilt_y, 2));
beam_tilt /= 1000.0f;
beam_tilt_x /= 1000.0f;
beam_tilt_y /= 1000.0f;
beam_tilt_z_X_component = tanf(beam_tilt) * cosf(beam_tilt_azimuth);
beam_tilt_z_Y_component = tanf(beam_tilt) * sinf(beam_tilt_azimuth);
}
// TODO fix periodicity in Z on slab
int iTilt;
int number_of_images;
float max_tilt = 0;
std::vector<float> tilt_psi;
std::vector<float> tilt_theta;
std::vector<float> tilt_phi;
if ( this->make_tilt_series ) {
number_of_images = n_tilt_angles;
if ( this->use_existing_params ) {
wxPrintf("\n\nUsing an existing parameter file only supported on a particle stack\n\n");
exit(-1);
}
// max_tilt = 60.0f;
// float tilt_range = max_tilt;
// float tilt_inc = 3.0f;
max_tilt = 0.0f;
for ( int iTilt = 0; iTilt < number_of_images; iTilt++ ) {
if ( fabsf(SET_TILT_ANGLES[iTilt]) > fabsf(max_tilt) ) {
max_tilt = SET_TILT_ANGLES[iTilt];
}
}
// float tilt_range = max_tilt;
float tilt_inc = 1.4;
for ( int iTilt = 0; iTilt < number_of_images; iTilt++ ) {
// to create a conical tilt (1.8*iTilt)+
tilt_psi.push_back(tilt_axis + this->stdErr * my_rand.GetNormalRandomSTD(0.0f, in_plane_sigma));
; // *(2*PI);
// tilt_theta[iTilt] = -((tilt_range - (float)iTilt*tilt_inc) + this->stdErr * my_rand.GetExponentialRandomSTD(0.0f,tilt_angle_sigma);
tilt_theta.push_back(SET_TILT_ANGLES[iTilt]);
if ( DO_PRINT ) {
wxPrintf("%f\n", SET_TILT_ANGLES[iTilt]);
}
tilt_phi.push_back(0.0f);
shift_x.push_back(this->stdErr * my_rand.GetUniformRandomSTD(-8.f, 8.f));
shift_y.push_back(this->stdErr * my_rand.GetUniformRandomSTD(-8.f, 8.f));
shift_z.push_back(0.0f);
mag_diff.push_back(1.0f + this->stdErr * my_rand.GetNormalRandomSTD(0.0f, magnification_sigma));
}
}
else if ( this->make_particle_stack > 0 ) {
// Set this up to get a restriction on angles due to symmetry.
ParameterMap my_dummy_map;
EulerSearch my_angle_restrictions;
my_angle_restrictions.InitGrid(wanted_symmetry, 5., 0., 0., 360, 1.5, 0., wanted_pixel_size * 2, my_dummy_map, 5);
float r_phi_max = my_angle_restrictions.phi_max;
float r_theta_max = (my_angle_restrictions.test_mirror) ? 2.f * my_angle_restrictions.theta_max : my_angle_restrictions.theta_max;
max_tilt = 0.0f;
number_of_images = this->make_particle_stack;
std::normal_distribution<float> normal_dist(0.0f, 1.0f);
tilt_psi.reserve(number_of_images);