AUTHORS

-This program is developed in the group of Sjors H.W. Scheres at the MRC Laboratory of Molecular Biology.
+This program is developed in the group of Sjors H.W. Scheres at the MRC Laboratory of Molecular Biology,
+
+- Sjors H.W. Scheres
+- Shaoda He
+- Takanori Nakane
+- Jasenko Zivanov
+- Liyi Dong
+
+and Erik Lindahl at Stockholm University.
+
+- Erik Lindahl
+- Björn O. Forsberg
+- Dari Kimanius
 
 However, it does also contain pieces of code from the following packages:
 XMIPP: http:/xmipp.cnb.csic.es

CMakeLists.txt

@@ -184,7 +184,7 @@ if(GUI)
 				message(STATUS "X11 and FLTK were found")
 				message(STATUS "FLTK_LIBRARIES: ${FLTK_LIBRARIES}")
 		    else()
-		        message(STATUS "FLTK was NOT found")
+				message(STATUS "No FLTK installation was found")
 		    endif() 
         endif(NOT FORCE_OWN_FLTK)
         
@@ -210,6 +210,7 @@ if(NOT FORCE_OWN_FFTW)
 endif(NOT FORCE_OWN_FFTW)
 	
 if(NOT FFTW_FOUND)
+	message(STATUS "No FFTW installation was found")
 	include(${CMAKE_SOURCE_DIR}/cmake/BuildFFTW.cmake)
 endif(NOT FFTW_FOUND)
 

README.md

-
-
-#RELION
+RELION
+======
 
 
 RELION (for REgularised LIkelihood OptimisatioN) is a stand-alone computer
@@ -22,10 +21,10 @@ More extensive options and configurations are available
 but the outlines to clone and install relion for typical use are made easy
 through [cmake](https://en.wikipedia.org/wiki/CMake).
 
-On ubuntu machines, installing cmake is as easy as 
+On ubuntu machines, installing cmake, the compiler, and additional dependencies (mpi, fftw) is as easy as:
 
 ```
-sudo apt install cmake
+sudo apt install cmake build-essential mpi-default-bin mpi-default-dev libfftw3-dev
 ```
 
 On other systems it is typically just as easy, you simply have to modify "apt" to
@@ -46,9 +45,13 @@ cmake -DCMAKE_INSTALL_PREFIX=/where/to/install/ ..
 make -j4
 make install
 ```
+(NOTES: "/where/to/install/.." is typically "/usr/local/relion".
+ Make sure you create this directory beforehand.
+ Installing to that location requires sudo, so in this case be sure to use
+ "sudo make install" instead of "make install" in the last step.)
 
 These steps will download the source-code, create a build-directory,
-then configure and build relion, and lastely install it to be generally
+then configure and build relion, and lastly install it to be generally
 available on the system.
 
 
@@ -59,9 +62,10 @@ RELION is intermittently updated, with both minor and major features.
 To update an existing installation, simply use the following commands
 
 ```
-cd relion/build
+cd relion
 git pull
+cd build
 make -j4
-make
+make install    # (or "sudo make install")
 
 ```

cmake/BuildFFTW.cmake

-message(STATUS "-------------------------------------------------")
-message(STATUS "------- WILL USE LOCALLY BUILT FFTW LIBS --------")
-message(STATUS "-------------------------------------------------")
-
 set(FFTW_EXTERNAL_PATH "${CMAKE_SOURCE_DIR}/external/fftw")
 
 if (NOT DEFINED TARGET_X86)
@@ -32,9 +28,16 @@ find_library(FFTW_LIBRARIES NAMES ${libfft}       PATHS ${FFTW_EXTERNAL_PATH}/li
 
 if(FFTW_INCLUDES AND FFTW_LIBRARIES)
     set(FFTW_FOUND TRUE)
-    message( STATUS "Found previously built external (non-system) FFTW library")
+	message(STATUS "Found previously built non-system FFTW libraries that will be used.")
 else()
     set(FFTW_FOUND FALSE)
+	message(STATUS "--------------------------------------------------------")
+	message(STATUS "-------- NO EXISTING FFTW LIBRARIES WHERE FOUND. -------")
+	message(STATUS "-------------- FFTW WILL BE DOWNLOADED AND -------------")
+	message(STATUS "--------------- BUILT DURING COMPILE-TIME. -------------")
+	message(STATUS "--------------------------------------------------------")
+	message(STATUS "---- A WORKING INTERNET CONNECTION WILL BE REQUIRED. ---")
+	message(STATUS "--------------------------------------------------------")
 endif()
 
 ## ----------------------------------------------------------------- NEW EXT LIBS? --

cmake/BuildFLTK.cmake

-message(STATUS "-------------------------------------------------")
-message(STATUS "------- WILL USE LOCALLY BUILT FLTK LIBS --------")
-message(STATUS "-------------------------------------------------")
+
 
 set(FLTK_EXTERNAL_PATH "${CMAKE_SOURCE_DIR}/external/fltk")
 set(ext_conf_flags_fltk --enable-shared --prefix=${FLTK_EXTERNAL_PATH})
@@ -13,10 +11,18 @@ find_path(FLTK_INCLUDES     NAMES FL/Fl.H   PATHS  "${FLTK_EXTERNAL_PATH}/includ
 
 if(FLTK_INCLUDE_DIR AND FLTK_LIBRARIES)
     set(FLTK_FOUND TRUE)
-    message( STATUS "Found previously built external (non-system) FLTK library")
+	message(STATUS "Found previously built non-system FLTK libraries that will be used.")
 else()
     set(FLTK_FOUND FALSE)
+	message(STATUS "--------------------------------------------------------")
+	message(STATUS "-------- NO EXISTING FLTK LIBRARIES WHERE FOUND. -------")
+	message(STATUS "-------------- FLTK WILL BE DOWNLOADED AND -------------")
+	message(STATUS "--------------- BUILT DURING COMPILE-TIME. -------------")
+	message(STATUS "--------------------------------------------------------")
+	message(STATUS "---- A WORKING INTERNET CONNECTION WILL BE REQUIRED. ---")
+	message(STATUS "--------------------------------------------------------")
 endif()
+
 ## ----------------------------------------------------------------- NEW EXT LIBS? --  
  
 if(NOT FLTK_FOUND)

debian/changelog

-relion (2.1.b1-1) UNRELEASED; urgency=medium
+relion (2.1-1) UNRELEASED; urgency=medium
 
   * Team upload.
   * Add watch file

debian/watch

 version=4
 
-https://github.com/3dem/relion/releases .*/archive/@ANY_VERSION@@ARCHIVE_EXT@
+https://github.com/3dem/relion/releases .*/archive/(\d[\d.]+)@ARCHIVE_EXT@

src/apps/CMakeLists.txt

@@ -8,31 +8,41 @@ find_path (X11_INCLUDES Xdbe.h)
 message(STATUS "CMAKE_BINARY_DIR:" ${CMAKE_BINARY_DIR})
 
 file(GLOB REL_SRC "${CMAKE_SOURCE_DIR}/src/*.cpp" "${CMAKE_SOURCE_DIR}/src/*.c" "${CMAKE_SOURCE_DIR}/src/gpu_utils/*.cpp")
+file(GLOB REL_GUI_SRC "${CMAKE_SOURCE_DIR}/src/manualpicker.cpp" "${CMAKE_SOURCE_DIR}/src/gui_*.cpp" "${CMAKE_SOURCE_DIR}/src/displayer.cpp")
+# Remove GUI files from relion_lib
+foreach(GUI_SRC_FILE ${REL_GUI_SRC})
+	list(REMOVE_ITEM REL_SRC "${GUI_SRC_FILE}")
+endforeach()
 file(GLOB REL_SRC_H "${CMAKE_SOURCE_DIR}/src/*.h" "${CMAKE_SOURCE_DIR}/src/gpu_utils/*.h")
 file(GLOB REL_HP "${CMAKE_SOURCE_DIR}/src/Healpix_2.15a/*.cc")
 
 file(GLOB RELION_TARGETS "${CMAKE_SOURCE_DIR}/src/apps/*.cpp")
 
+set(GUI_TARGETS display maingui manualpick pipeliner)
+
 #--Remove apps using X11 if no GUI--
 if(NOT GUI)
-    list(REMOVE_ITEM REL_SRC "${CMAKE_SOURCE_DIR}/src/manualpicker.cpp")
-    list(REMOVE_ITEM REL_SRC "${CMAKE_SOURCE_DIR}/src/gui_entries.cpp")
-    list(REMOVE_ITEM REL_SRC "${CMAKE_SOURCE_DIR}/src/gui_jobwindow.cpp")
-    list(REMOVE_ITEM REL_SRC "${CMAKE_SOURCE_DIR}/src/gui_mainwindow.cpp")
-    list(REMOVE_ITEM REL_SRC "${CMAKE_SOURCE_DIR}/src/displayer.cpp")
-    list(REMOVE_ITEM RELION_TARGETS "${CMAKE_SOURCE_DIR}/src/apps/display.cpp")
-    list(REMOVE_ITEM RELION_TARGETS "${CMAKE_SOURCE_DIR}/src/apps/maingui.cpp")
-    list(REMOVE_ITEM RELION_TARGETS "${CMAKE_SOURCE_DIR}/src/apps/manualpick.cpp")
-    list(REMOVE_ITEM RELION_TARGETS "${CMAKE_SOURCE_DIR}/src/apps/pipeliner.cpp")
+    foreach(TARGET ${GUI_TARGETS})
+        list(REMOVE_ITEM RELION_TARGETS "${CMAKE_SOURCE_DIR}/src/apps/${TARGET}.cpp")
+    endforeach()
 endif(NOT GUI)
 
 
 # relion_lib is STATIC or SHARED type based on BUILD_SHARED_LIBS=ON/OFF
+# relion_lib only contains non-X11 parts
+# relion_gui_lib is where the X11 code is placed
 if(BUILD_SHARED_LIBS)
     add_library(relion_lib SHARED ${REL_SRC} ${REL_SRC_H} ${REL_HP})
     install(TARGETS relion_lib LIBRARY DESTINATION lib)
+    if(GUI)
+        add_library(relion_gui_lib SHARED ${REL_GUI_SRC} ${REL_SRC_H} ${REL_HP})
+        install(TARGETS relion_gui_lib LIBRARY DESTINATION lib)
+    endif(GUI)
 else()
     add_library(relion_lib STATIC ${REL_SRC} ${REL_SRC_H} ${REL_HP})
+    if(GUI)
+        add_library(relion_gui_lib STATIC ${REL_GUI_SRC} ${REL_SRC_H} ${REL_HP})
+    endif(GUI)
 endif()
 
 target_link_libraries(relion_lib ${FFTW_LIBRARIES})
@@ -41,9 +51,10 @@ if(BUILD_OWN_FFTW)
 endif()
 
 if(GUI)
-    target_link_libraries(relion_lib ${FLTK_LIBRARIES})
+    include_directories("${FLTK_INCLUDE_DIR}")
+    target_link_libraries(relion_gui_lib relion_lib ${FLTK_LIBRARIES})
     if(BUILD_OWN_FLTK)
-        add_dependencies(relion_lib OWN_FLTK)
+        add_dependencies(relion_gui_lib OWN_FLTK)
     endif()
 endif(GUI)
 
@@ -101,10 +112,16 @@ foreach (_target ${RELION_TARGETS})
 	    target_link_libraries(${_target} relion_gpu_util)
 	endif(CUDA_FOUND)
 	
-	if(GUI)
-        include_directories("${FLTK_INCLUDE_DIR}")
-		target_link_libraries(${_target} ${FLTK_LIBRARIES} ${X11})
-	endif(GUI)
+#	if(${_target} IN_LIST GUI_TARGETS)
+#        add_dependencies(${_target} relion_gui_lib)
+#		target_link_libraries(${_target} relion_gui_lib ${FLTK_LIBRARIES} ${X11})
+#	endif()
+	
+	list(FIND GUI_TARGETS ${_target} IS_GUI_TARGET)
+	if(NOT ${IS_GUI_TARGET} LESS 0)
+        add_dependencies(${_target} relion_gui_lib)
+        target_link_libraries(${_target} relion_gui_lib ${FLTK_LIBRARIES} ${X11})
+    endif()
     
 	message(STATUS "added ${_target}...")
 	
@@ -120,3 +137,4 @@ if(GUI)
 endif(GUI)
 
 
+

src/apps/helix_toolbox.cpp

@@ -101,6 +101,9 @@ public:
 	// Rotational symmetry - Cn
 	int sym_Cn;
 
+	// Number of filaments in a helix with seam (>= 2)
+	int nr_filaments_helix_with_seam;
+
 	// Helical rise and its local searches
 	RFLOAT rise_A, rise_min_A, rise_max_A, rise_inistep_A;
 
@@ -321,6 +324,7 @@ public:
 		rise_inistep_A = textToFloat(parser.getOption("--rise_inistep", "Initial step of helical rise search (in Angstroms)", "-1"));
 		rise_min_A = textToFloat(parser.getOption("--rise_min", "Minimum helical rise (in Angstroms)", "-1"));
 		rise_max_A = textToFloat(parser.getOption("--rise_max", "Maximum helical rise (in Angstroms)", "-1"));
+		nr_filaments_helix_with_seam = textToInteger(parser.getOption("--seam_nr_filaments", "Number of filaments in a helix with seam (>= 2)", "-1"));
 		do_helical_symmetry_local_refinement = parser.checkOption("--search_sym", "Perform local searches of helical symmetry in 3D reconstruction?");
 		do_cut_into_segments = parser.checkOption("--segments", "Cut helical tubes into segments?");
 		sigma_offset = textToFloat(parser.getOption("--sigma_offset", "Sigma of translational offsets (in pixels)", "5."));
@@ -347,6 +351,10 @@ public:
 			SWAP(rise_min_A, rise_max_A, tmp_RFLOAT);
 		if (twist_min_deg > twist_max_deg)
 			SWAP(twist_min_deg, twist_max_deg, tmp_RFLOAT);
+
+		// Check for errors in the command-line option
+		if (parser.checkForErrors())
+			REPORT_ERROR("Errors encountered on the command line (see above), exiting...");
 	};
 
 	void clear()
@@ -627,7 +635,7 @@ public:
 			{
 				displayEmptyLine();
 				std::cout << " Create a helical 3D reference of spheres" << std::endl;
-				std::cout << "  USAGE: --simulate_helix --o ref.mrc --subunit_diameter 30 --cyl_outer_diameter 200 --angpix 1.126 --rise 1.408 --twist 22.03 --boxdim 300 (--sym_Cn 1) (--polar --topbottom_ratio 0.5 --cyl_inner_diameter 20) (--sphere_percentage 0.9 --width 5)" << std::endl;
+				std::cout << "  USAGE: --simulate_helix --o ref.mrc --subunit_diameter 30 --cyl_outer_diameter 200 --angpix 1.126 --rise 1.408 --twist 22.03 --boxdim 300 (--sym_Cn 1) (--polar --topbottom_ratio 0.5 --cyl_inner_diameter 20) (--sphere_percentage 0.9 --width 5) (--seam_nr_filaments 13)" << std::endl;
 				displayEmptyLine();
 				return;
 			}
@@ -650,7 +658,8 @@ public:
 					subunit_diameter_A,
 					(do_polar_reference) ? (cyl_inner_diameter_A) : (subunit_diameter_A),
 					(do_polar_reference) ? (topbottom_ratio) : (1.),
-					sym_Cn);
+					sym_Cn,
+					nr_filaments_helix_with_seam);
 			applySoftSphericalMask(
 					img(),
 					(RFLOAT(boxdim) * sphere_percentage),

src/apps/image_handler.cpp

@@ -30,7 +30,7 @@ class image_handler_parameters
 	public:
    	FileName fn_in, fn_out, fn_sel, fn_img, fn_sym, fn_sub, fn_mult, fn_div, fn_add, fn_subtract, fn_fsc, fn_adjust_power, fn_correct_ampl, fn_fourfilter;
 	int bin_avg, avg_first, avg_last, edge_x0, edge_xF, edge_y0, edge_yF, filter_edge_width, new_box, minr_ampl_corr;
-    bool do_add_edge, do_flipXY, do_flipmXY, do_flipZ, do_flipX, do_flipY, do_shiftCOM, do_stats, do_avg_ampl, do_average, do_remove_nan;
+    bool do_add_edge, do_flipXY, do_flipmXY, do_flipZ, do_flipX, do_flipY, do_shiftCOM, do_stats, do_avg_ampl, do_avg_ampl2, do_avg_ampl2_ali, do_average, do_remove_nan;
 	RFLOAT multiply_constant, divide_constant, add_constant, subtract_constant, threshold_above, threshold_below, angpix, new_angpix, lowpass, highpass, bfactor, shift_x, shift_y, shift_z, replace_nan;
    	int verb;
 	// I/O Parser
@@ -94,6 +94,8 @@ class image_handler_parameters
 	    shift_y = textToFloat(parser.getOption("--shift_y", "Shift images this many pixels in the Y-direction", "0."));
 	    shift_z = textToFloat(parser.getOption("--shift_z", "Shift images this many pixels in the Z-direction", "0."));
 	    do_avg_ampl = parser.checkOption("--avg_ampl", "Calculate average amplitude spectrum for all images?");
+	    do_avg_ampl2 = parser.checkOption("--avg_ampl2", "Calculate average amplitude spectrum for all images?");
+	    do_avg_ampl2_ali = parser.checkOption("--avg_ampl2_ali", "Calculate average amplitude spectrum for all aligned images?");
 	    do_average = parser.checkOption("--average", "Calculate average of all images (without alignment)");
 	    fn_correct_ampl = parser.getOption("--correct_avg_ampl", "Correct all images with this average amplitude spectrum", "");
 	    minr_ampl_corr = textToInteger(parser.getOption("--minr_ampl_corr", "Minimum radius (in Fourier pixels) to apply average amplitudes", "0"));
@@ -423,7 +425,7 @@ class image_handler_parameters
 		// Thresholding (can be done after any other operation)
 		if (fabs(threshold_above - 999.) > 0.)
 		{
-			FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(Iin())
+			FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(Iout())
 			{
 				if (DIRECT_A3D_ELEM(Iout(), k, i, j) > threshold_above)
 					DIRECT_A3D_ELEM(Iout(), k, i, j) = threshold_above;
@@ -431,7 +433,7 @@ class image_handler_parameters
 		}
 		if (fabs(threshold_below + 999.) > 0.)
 		{
-			FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(Iin())
+			FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY3D(Iout())
 			{
 				if (DIRECT_A3D_ELEM(Iout(), k, i, j) < threshold_below)
 					DIRECT_A3D_ELEM(Iout(), k, i, j) = threshold_below;
@@ -584,7 +586,7 @@ class image_handler_parameters
 					if (XSIZE(Iop()) != xdim || YSIZE(Iop()) != ydim || ZSIZE(Iop()) != zdim)
 						REPORT_ERROR("Error: operate-image is not of the correct size");
 
-				if (do_avg_ampl)
+				if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali)
 				{
 					avg_ampl.initZeros(zdim, ydim, xdim/2+1);
 				}
@@ -603,16 +605,44 @@ class image_handler_parameters
 				Iin().computeStats(avg, stddev, minval, maxval);
 				std::cout << fn_img << " : (x,y,z,n)= " << XSIZE(Iin()) << " x "<< YSIZE(Iin()) << " x "<< ZSIZE(Iin()) << " x "<< NSIZE(Iin()) << " ; avg= " << avg << " stddev= " << stddev << " minval= " <<minval << " maxval= " << maxval << std::endl;
 			}
-			else if (do_avg_ampl)
+			else if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali)
 			{
 				Iin.read(fn_img);
+
+				if (do_avg_ampl2_ali)
+				{
+					RFLOAT xoff = 0.;
+					RFLOAT yoff = 0.;
+					RFLOAT psi = 0.;
+					MD.getValue(EMDL_ORIENT_ORIGIN_X, xoff);
+					MD.getValue(EMDL_ORIENT_ORIGIN_Y, yoff);
+					MD.getValue(EMDL_ORIENT_PSI, psi);
+					// Apply the actual transformation
+					Matrix2D<RFLOAT> A;
+					rotation2DMatrix(psi, A);
+				    MAT_ELEM(A,0, 2) = xoff;
+				    MAT_ELEM(A,1, 2) = yoff;
+				    selfApplyGeometry(Iin(), A, IS_NOT_INV, DONT_WRAP);
+				}
+
 				MultidimArray<Complex> FT;
 				transformer.FourierTransform(Iin(), FT);
+
+				if (do_avg_ampl)
+				{
 					FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FT)
 					{
 						DIRECT_MULTIDIM_ELEM(avg_ampl, n) +=  abs(DIRECT_MULTIDIM_ELEM(FT, n));
 					}
 				}
+				else if (do_avg_ampl2 || do_avg_ampl2_ali)
+				{
+					FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(FT)
+					{
+						DIRECT_MULTIDIM_ELEM(avg_ampl, n) +=  norm(DIRECT_MULTIDIM_ELEM(FT, n));
+					}
+				}
+			}
 			else if (do_average)
 			{
 				Iin.read(fn_img);
@@ -668,6 +698,10 @@ class image_handler_parameters
 				FileName my_fn_out;
 				if (fn_out == "")
 					my_fn_out = fn_img;
+				else if(fn_out.getExtension() == "mrcs" && !fn_out.contains("@"))
+				{
+					my_fn_out.compose(current_object,fn_out);
+				}
 				else
 				{
 					if (fn_in.getExtension() == "star" || (fn_in.getExtension() == "mrcs" && !fn_in.contains("@")))
@@ -690,7 +724,7 @@ class image_handler_parameters
 		}
 
 
-		if (do_avg_ampl || do_average)
+		if (do_avg_ampl || do_avg_ampl2 || do_avg_ampl2_ali || do_average)
 		{
 			avg_ampl /= (RFLOAT)i_img;
 			Iout() = avg_ampl;

src/apps/mask_create.cpp

@@ -62,12 +62,13 @@ class mask_create_parameters
 	    angpix = textToFloat(parser.getOption("--angpix", "Pixel size (in Angstroms) for the lowpass filter", "1"));
 	    helical_z_percentage = textToFloat(parser.getOption("--z_percentage", "This box length along the center of Z axis contains good information of the helix", "0.3"));
 
-	    if (fn_thr=="" && fn_apply_in == "")
-	    	REPORT_ERROR("Either provide --i to apply a mask, OR --create_from to create a new mask");
-
 	    // Check for errors in the command-line option
     	if (parser.checkForErrors())
     		REPORT_ERROR("Errors encountered on the command line, exiting...");
+
+    	if (fn_thr=="" && fn_apply_in == "")
+	    	REPORT_ERROR("Either provide --i to apply a mask, OR --create_from to create a new mask");
+
 	}
 
 

src/apps/particle_symmetry_expand.cpp

@@ -65,6 +65,10 @@ public:
         frac_sampling = textToFloat(parser.getOption("--frac_sampling", "Number of samplings in between a single asymmetrical unit", "1"));
         frac_range = textToFloat(parser.getOption("--frac_range", "Range of the rise [-0.5, 0.5> to be sampled", "0.5"));
 
+       	// Check for errors in the command-line option
+    	if (parser.checkForErrors())
+    		REPORT_ERROR("Errors encountered on the command line (see above), exiting...");
+
         if (do_helix)
         {
         	if (fn_sym != "C1")
@@ -75,10 +79,6 @@ public:
 
         }
 
-       	// Check for errors in the command-line option
-    	if (parser.checkForErrors())
-    		REPORT_ERROR("Errors encountered on the command line (see above), exiting...");
-
 	}
 
 	void run()

src/apps/project.cpp

@@ -37,9 +37,9 @@ class project_parameters
 public:
 
 	FileName fn_map, fn_ang, fn_out, fn_img, fn_model, fn_sym, fn_mask;
-	RFLOAT rot, tilt, psi, xoff, yoff, zoff, angpix, maxres, stddev_white_noise, particle_diameter, ana_prob_range, ana_prob_step, sigma_offset;
+	RFLOAT rot, tilt, psi, xoff, yoff, zoff, angpix, maxres, stddev_white_noise, particle_diameter, ana_prob_range, ana_prob_step;
 	int padding_factor;
-	int r_max, r_min_nn, interpolator, nr_uniform;
+	int r_max, r_min_nn, interpolator;
     bool do_only_one, do_ctf, do_ctf2, ctf_phase_flipped, do_ctf_intact_1st_peak, do_timing, do_add_noise, do_subtract_exp, do_ignore_particle_name, do_3d_rot;
 	// I/O Parser
 	IOParser parser;
@@ -64,8 +64,6 @@ public:
        	angpix = textToFloat(parser.getOption("--angpix", "Pixel size (in Angstroms)", "1"));
 		fn_mask = parser.getOption("--mask", "Mask that will be applied to the input map prior to making projections", "");
        	fn_ang = parser.getOption("--ang", "STAR file with orientations for multiple projections (if None, assume single projection)","None");
-       	nr_uniform = textToInteger(parser.getOption("--nr_uniform", " OR get this many samples from a uniform angular distribution", "-1"));
-       	sigma_offset = textToFloat(parser.getOption("--sigma_offset", "Apply Gaussian errors with this stddev to the XY-offsets", "0"));
        	rot = textToFloat(parser.getOption("--rot", "First Euler angle (for a single projection)", "0"));
        	tilt = textToFloat(parser.getOption("--tilt", "Second Euler angle (for a single projection)", "0"));
        	psi = textToFloat(parser.getOption("--psi", "Third Euler angle (for a single projection)", "0"));
@@ -77,7 +75,7 @@ public:
        	fn_model = parser.getOption("--model_noise", "Model STAR file with power spectra for coloured Gaussian noise", "");
        	do_subtract_exp = parser.checkOption("--subtract_exp", "Subtract projections from experimental images (in --ang)");
        	do_ignore_particle_name = parser.checkOption("--ignore_particle_name", "Ignore the rlnParticleName column (in --ang)");
-       	do_only_one = (fn_ang == "None" && nr_uniform < 0);
+       	do_only_one = (fn_ang == "None");
        	do_3d_rot = parser.checkOption("--3d_rot", "Perform 3D rotations instead of projection into 2D images");
 
        	maxres = textToFloat(parser.getOption("--maxres", "Maximum resolution (in Angstrom) to consider in Fourier space (default Nyquist)", "-1"));
@@ -123,35 +121,7 @@ public:
     		DIRECT_MULTIDIM_ELEM(vol(), n) *= DIRECT_MULTIDIM_ELEM(msk(), n);
     	}
 
-    	if (nr_uniform > 0)
-    	{
-
-    		std::cout << " Generating " << nr_uniform << " projections taken randomly from a uniform angular distribution ..." << std::endl;
-    		MDang.clear();
-    		randomize_random_generator();
-    		for (long int i = 0; i < nr_uniform; i++)
-    		{
-				RFLOAT rot, tilt, psi, xoff, yoff;
-				rot = rnd_unif() * 360.;
-				bool ok_tilt = false;
-				while (!ok_tilt)
-				{
-					tilt = rnd_unif() * 180.;
-					if (rnd_unif() < fabs(SIND(tilt)))
-						ok_tilt = true;
-				}
-				psi = rnd_unif() * 360.;
-				xoff = rnd_gaus(0., sigma_offset);
-				yoff = rnd_gaus(0., sigma_offset);
-				MDang.addObject();
-				MDang.setValue(EMDL_ORIENT_ROT, rot);
-				MDang.setValue(EMDL_ORIENT_TILT, tilt);
-				MDang.setValue(EMDL_ORIENT_PSI, psi);
-				MDang.setValue(EMDL_ORIENT_ORIGIN_X, xoff);
-				MDang.setValue(EMDL_ORIENT_ORIGIN_Y, yoff);
-    		}
-    	}
-    	else if (!do_only_one)
+    	if (!do_only_one)
     	{
     		std::cout << " Reading STAR file with all angles " << fn_ang << std::endl;
     		MDang.read(fn_ang);
@@ -449,7 +419,7 @@ int main(int argc, char *argv[])
 
     catch (RelionError XE)
     {
-        prm.usage();
+        //prm.usage();
         std::cerr << XE;
         exit(1);
     }

src/apps/reconstruct.cpp

@@ -33,7 +33,8 @@ class reconstruct_parameters
 	MetaDataTable DF;
 	int r_max, r_min_nn, blob_order, ref_dim, interpolator, iter, nr_threads, debug_ori_size, debug_size, ctf_dim, nr_helical_asu;
 	RFLOAT blob_radius, blob_alpha, angular_error, shift_error, angpix, maxres, beamtilt_x, beamtilt_y, helical_rise, helical_twist;
-	bool do_ctf, ctf_phase_flipped, only_flip_phases, intact_ctf_first_peak, do_fom_weighting, do_3d_rot, do_reconstruct_ctf, do_beamtilt, skip_gridding, do_reconstruct_ctf2;
+	bool do_ctf, ctf_phase_flipped, only_flip_phases, intact_ctf_first_peak, do_fom_weighting, do_3d_rot, do_reconstruct_ctf, do_beamtilt;
+	bool skip_gridding, do_reconstruct_ctf2, do_reconstruct_meas;
 	float padding_factor;
 	// I/O Parser
 	IOParser parser;
@@ -49,11 +50,6 @@ class reconstruct_parameters
 		parser.setCommandLine(argc, argv);
 
 		int general_section = parser.addSection("General options");
-
-		fn_debug = parser.getOption("--debug", "Rootname for debug reconstruction files", "");
-		debug_ori_size =  textToInteger(parser.getOption("--debug_ori_size", "Rootname for debug reconstruction files", "1"));
-		debug_size =  textToInteger(parser.getOption("--debug_size", "Rootname for debug reconstruction files", "1"));
-
 		fn_sel = parser.getOption("--i", "Input STAR file with the projection images and their orientations", "");
 	    fn_out = parser.getOption("--o", "Name for output reconstruction","relion.mrc");
 	    fn_sym = parser.getOption("--sym", "Symmetry group", "c1");
@@ -94,10 +90,15 @@ class reconstruct_parameters
     	do_3d_rot = parser.checkOption("--3d_rot", "Perform 3D rotations instead of backprojections from 2D images");
     	ctf_dim  = textToInteger(parser.getOption("--reconstruct_ctf", "Perform a 3D reconstruction from 2D CTF-images, with the given size in pixels", "-1"));
     	do_reconstruct_ctf2 = parser.checkOption("--ctf2", "Reconstruct CTF^2 and then take the sqrt of that");
+    	do_reconstruct_meas = parser.checkOption("--measured", "Fill Hermitian half of the CTF reconstruction with how often each voxel was measured.");
     	skip_gridding = parser.checkOption("--skip_gridding", "Skip gridding part of the reconstruction");
     	do_reconstruct_ctf = (ctf_dim > 0);
     	if (do_reconstruct_ctf)
     		do_ctf = false;
+		fn_debug = parser.getOption("--debug", "Rootname for debug reconstruction files", "");
+		debug_ori_size =  textToInteger(parser.getOption("--debug_ori_size", "Rootname for debug reconstruction files", "1"));
+		debug_size =  textToInteger(parser.getOption("--debug_size", "Rootname for debug reconstruction files", "1"));
+
 
     	// Hidden
        	r_min_nn = textToInteger(getParameter(argc, argv, "--r_min_nn", "10"));
@@ -243,6 +244,8 @@ class reconstruct_parameters
     	}
 
    		BackProjector backprojector(mysize, ref_dim, fn_sym, interpolator, padding_factor, r_min_nn, blob_order, blob_radius, blob_alpha, data_dim, skip_gridding);
+   		// This one is only needed for do_reconstruct_meas, but non-empty constructor does not exist...
+   		BackProjector backprojector2(mysize, ref_dim, fn_sym, interpolator, padding_factor, r_min_nn, blob_order, blob_radius, blob_alpha, data_dim, skip_gridding);
 
    		if (maxres < 0.)
    			r_max = -1;
@@ -259,7 +262,7 @@ class reconstruct_parameters
    		}
 
 		// Check for beam-tilt parameters in the input star file
-   		if (do_beamtilt && ( DF.containsLabel(EMDL_IMAGE_BEAMTILT_X) || DF.containsLabel(EMDL_IMAGE_BEAMTILT_Y) ) )
+   		if (do_beamtilt || ( DF.containsLabel(EMDL_IMAGE_BEAMTILT_X) || DF.containsLabel(EMDL_IMAGE_BEAMTILT_Y) ) )
    				std::cout << " + Using the beamtilt parameters in the input STAR file" << std::endl;
 
 		std::cerr << "Back-projecting all images ..." << std::endl;
@@ -433,6 +436,17 @@ class reconstruct_parameters
 
 
 				backprojector.set2DFourierTransform(F2D, A3D, IS_NOT_INV, &Fctf);
+
+
+				if (do_reconstruct_meas)
+				{
+					// Reconstruct from all-one F2Ds
+					FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(F2D)
+					{
+						DIRECT_MULTIDIM_ELEM(F2D, n) = 1.;
+					}
+					backprojector.set2DFourierTransform(F2D, A3D, IS_NOT_INV, &Fctf);
+				}
 			}
 
 
@@ -502,6 +516,57 @@ class reconstruct_parameters
 				}
    			}
 
+
+   			if (do_reconstruct_meas)
+   			{
+   				backprojector2.reconstruct(vol(), iter, do_map, 1., dummy, dummy, dummy, dummy, fsc, 1., do_use_fsc, true, nr_threads, -1);
+   	   			F2D.clear();
+   	   			transformer.FourierTransform(vol(), F2D);
+
+   	   			// CenterOriginFFT: Set the center of the FFT in the FFTW origin
+   				Matrix1D<RFLOAT> shift(3);
+   				XX(shift)=-(RFLOAT)(int)(ctf_dim / 2);
+   				YY(shift)=-(RFLOAT)(int)(ctf_dim / 2);
+   				ZZ(shift)=-(RFLOAT)(int)(ctf_dim / 2);
+   				shiftImageInFourierTransform(F2D, F2D, (RFLOAT)ctf_dim, XX(shift), YY(shift), ZZ(shift));
+
+   				// Divide each value by the average in each resolution shell.
+   				// This is the multiplicative correction that will later be applied to sigma2_noise
+   				MultidimArray<RFLOAT> sum, count;
+   				sum.initZeros(ctf_dim / 2);
+   				count.initZeros(ctf_dim / 2);
+   				FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM(F2D)
+   	   			{
+					int ires = ROUND(sqrt(kp*kp + ip*ip + jp*jp));
+					DIRECT_A1D_ELEM(sum, ires) += DIRECT_A3D_ELEM(F2D, k, i, j).real;
+					DIRECT_A1D_ELEM(count, ires) += 1.;
+   	   			}
+   				FOR_ALL_DIRECT_ELEMENTS_IN_ARRAY1D(sum)
+   				{
+   					if (DIRECT_A1D_ELEM(count, i) > 0.)
+   						DIRECT_A1D_ELEM(sum, i) /= DIRECT_A1D_ELEM(count, i);
+   				}
+  				FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM(F2D)
+   	   			{
+					int ires = ROUND(sqrt(kp*kp + ip*ip + jp*jp));
+   					if (DIRECT_A1D_ELEM(count, ires) > 0.)
+   						DIRECT_A3D_ELEM(F2D, k, i, j) /=  DIRECT_A1D_ELEM(sum, ires);
+   	   			}
+
+   				FOR_ALL_ELEMENTS_IN_FFTW_TRANSFORM(F2D)
+   	   			{
+   	   				// Take care of kp==dim/2, as XmippOrigin lies just right off center of image...
+   	   				if ( kp > FINISHINGZ(vol()) || ip > FINISHINGY(vol()) || jp > FINISHINGX(vol()))
+   	   					continue;
+   	   				if (kp == 0 && ip == 0 && jp == 0)
+   	   					continue;
+   	   				//A3D_ELEM(vol(), kp, ip, jp)    = FFTW_ELEM(F2D, kp, ip, jp).real;
+   	   				//NOW only set one Hermitian half to the number of measured components!
+   	   				A3D_ELEM(vol(), -kp, -ip, -jp) = FFTW_ELEM(F2D, kp, ip, jp).real;
+   	   			}
+   	   			vol() *= (RFLOAT)ctf_dim;
+   			}
+
    		}
 
    		vol.write(fn_out);

src/apps/stack_create.cpp

@@ -131,7 +131,7 @@ class stack_create_parameters
 
 			// Resize the output image
 			std::cout << "Resizing the output stack to "<< ndim<<" images of size: "<<xdim<<"x"<<ydim<<"x"<<zdim << std::endl;
-			RFLOAT Gb = (ndim*zdim*ydim*xdim*sizeof(RFLOAT))/(1024.*1024.*1024.);
+			RFLOAT Gb = (RFLOAT)ndim * zdim * ydim * xdim * sizeof(RFLOAT) / (1024. * 1024. * 1024.);
 			std::cout << "This will require " << Gb << "Gb of memory...."<< std::endl;
 			Image<RFLOAT> out(xdim, ydim, zdim, ndim);
 

src/args.cpp

@@ -249,6 +249,11 @@ void IOParser::writeCommandLine(std::ostream &out)
 bool IOParser::checkForErrors(int verb)
 {
 
+	if(checkParameter(argc, argv, "--version"))
+	{
+		std::cout << "RELION version " << RELION_VERSION << std::endl;
+		exit(0);
+	}
     if(argc==1 || (argc==2 && checkParameter(argc, argv, "--continue")) || checkParameter(argc, argv, "--help") || checkParameter(argc, argv, "-h"))
 	{
 		writeUsage(std::cout);

src/backprojector.cpp

@@ -901,6 +901,8 @@ void BackProjector::reconstruct(MultidimArray<RFLOAT> &vol_out,
 			}
 			myfsc = XMIPP_MIN(0.999, myfsc);
 			RFLOAT myssnr = myfsc / (1. - myfsc);
+			// Sjors 29nov2017 try tau2_fudge for pulling harder on Refine3D runs...
+			myssnr *= tau2_fudge;
 			RFLOAT fsc_based_tau = myssnr * DIRECT_A1D_ELEM(sigma2, i);
 			DIRECT_A1D_ELEM(tau2, i) = fsc_based_tau;
 			// data_vs_prior is merely for reporting: it is not used for anything in the reconstruction

src/ctf.cpp

@@ -175,8 +175,6 @@ void CTF::initialise()
     K3 = sqrt(1-Q0*Q0);
 
     K4 = -Bfac / 4.;
-    // For now switch off B-factor weighting on the CTF!
-    //K4 = 0.;
 
     // Phase shift in radian
     K5 = DEG2RAD(phase_shift);

src/ctffind_runner.cpp

@@ -67,7 +67,7 @@ void CtffindRunner::read(int argc, char **argv, int rank)
 	phase_max  = textToFloat(parser.getOption("--phase_max", "Maximum phase shift (in degrees)", "180."));
 	phase_step = textToFloat(parser.getOption("--phase_step", "Step in phase shift (in degrees)", "10."));
 	nr_threads = textToInteger(parser.getOption("--j", "Number of threads (for CTFIND4 only)", "1"));
-	do_large_astigmatism = parser.checkOption("--large_astigmatism", "Do you expect large astigmatism?");
+	do_fast_search = parser.checkOption("--fast_search", "Disable \"Slower, more exhaustive search\" in CTFFIND4.1 (faster but less accurate)");
 
 	int gctf_section = parser.addSection("Gctf parameters");
 	do_use_gctf = parser.checkOption("--use_gctf", "Use Gctf instead of CTFFIND to estimate the CTF parameters");
@@ -585,8 +585,12 @@ void CtffindRunner::executeCtffind4(long int imic)
 	fh << step_defocus << std::endl;
 	// Do you know what astigmatism is present?
 	fh << "no" << std::endl;
-	// Do you expect very large astigmatism?
-	if (do_large_astigmatism)
+	// Slower, more exhaustive search?
+	// The default was "no" in CTFFIND 4.1.5, but turned out to be less accurate.
+	// The default was changed to "yes" in CTFFIND 4.1.8. 
+	// Ref: http://grigoriefflab.janelia.org/ctffind4
+	// So, we say "yes" regardless of the version unless "--fast_search" is specified.
+	if (!do_fast_search)
 		fh << "yes" << std::endl;
 	else
 		fh << "no" << std::endl;

src/ctffind_runner.h

@@ -135,8 +135,8 @@ public:
 	// Current working directory to make absolute-path symlinks
 	std::string currdir;
 
-	// Expect large astigmatism for CTFFIND4.1?
-	bool do_large_astigmatism;
+	// Disable "Slower, more exhaustive search?" in CTFFIND 4.1.5-
+	bool do_fast_search;
 
 	// Which GPU devices to use?
 	std::string gpu_ids;