Additions to kernel

kernel/executables/2dx_backproject.cpp

@@ -0,0 +1,189 @@
+/* 
+ * @license GNU Public License
+ * @author Nikhil Biyani ([email protected])
+ * 
+ */
+
+#include <iostream>
+#include <string.h>
+
+#include "2dx_toolkit.h"
+
+template<typename ValueType_>
+std::vector<ValueType_>
+matrix_multiply(const std::vector<ValueType_>& mat1, const std::vector<ValueType_>& mat2, int r1, int c1, int c2){
+
+    std::vector<ValueType_> mult(r1*c2, ValueType_());
+
+    for(int i = 0; i < r1; ++i) {
+        for(int j = 0; j < c2; ++j) {
+            for(int k = 0; k < c1; ++k) {
+                mult[i+j*r1] += mat1[i+k*r1] * mat2[k+j*c1];
+            }
+        }
+    }
+
+    return mult;
+
+}
+
+void add_index(const tdx::data::MillerIndex& index_in, const tdx::Complex& value, 
+        tdx::data::MillerToPeakMultiMap& current_refs, double psi, double theta, double phi) {
+
+    // Values for transformation
+    double cpsi = cos(psi);
+    double cthe = cos(theta);
+    double cphi = cos(phi);
+    double spsi = sin(psi);
+    double sthe = sin(theta);
+    double sphi = sin(phi);
+
+    //TRANFORMATION MATRICES
+    std::vector<double> A = {cpsi, -1*spsi, 0, spsi, cpsi, 0, 0, 0, 1};
+    std::vector<double> B = {cthe, 0, sthe, 0, 1, 0, -1*sthe, 0, cthe};
+    std::vector<double> C = {cphi, -1*sphi, 0, sphi, cphi, 0, 0, 0, 1};
+    std::vector<double> in = {index_in.h()*1.0, index_in.k()*1.0, index_in.l()*1.0};
+
+    std::vector<double> res = matrix_multiply(in, matrix_multiply(A, matrix_multiply(B, C, 3, 3, 3), 3, 3, 3), 1, 3, 3);
+
+    int h_min = floor(res.at(0));
+    int k_min = floor(res.at(1));
+    int l_min = floor(res.at(2));
+
+    int h_max = ceil(res.at(0));
+    int k_max = ceil(res.at(1));
+    int l_max = ceil(res.at(2));
+
+    for(int h = h_min; h <= h_max; ++h) {
+        for(int k = k_min; k<= k_max; ++k) {
+            for(int l = l_min; l<-l_max; ++l) {
+
+                double dist = sqrt((res[0]-h)*(res[0]-h) + (res[1]-k)*(res[1]-k) + (res[2]-l)*(res[2]-l));
+                double sinc = 1;
+                if(dist != 0) sinc = sin(M_PI*dist)/(M_PI*dist);
+                auto amp_in = value.amplitude();
+                auto phase_in = value.phase();
+
+                tdx::data::MillerIndex new_idx(h, k, l);
+
+                if(new_idx.h() < 0) {
+                    new_idx = new_idx.FriedelSpot();
+                    phase_in *= -1;
+                }
+
+                double real_in = amp_in*sinc*cos(phase_in);
+                double imag_in = amp_in*sinc*sin(phase_in);
+
+                tdx::Complex complex_in(real_in, imag_in);
+                tdx::data::PeakData value_out(complex_in, 1.0);
+
+                current_refs.insert(tdx::data::MillerToPeakPair(new_idx, value_out));
+            }
+        }
+    }
+}
+
+int main(int argc, char* argv[])
+{
+
+    args::Executable exe("Program to backproject particles to 3D volume.", ' ', "1.0" );
+
+    //Select required arguments
+    TCLAP::ValueArg<std::string> INPUT_VOLUME("", "input", "Input MRC stack file", true, "", "MRC FILE");
+    TCLAP::ValueArg<std::string> INPUT_PAR("", "par", "Input PAR file", true, "", "PAR FILE");
+    TCLAP::ValueArg<std::string> OUTPUT_VOLUME("", "output", "OUTPUT Volume file", true, "", "MRC FILE");
+    TCLAP::ValueArg<int> NUM_PARTICLES("", "particles", "Number of particles to consider", false, 0, "INT");
+    TCLAP::ValueArg<double> APIX("", "apix", "Pixel size (A/pixel)", false, 1.0, "FLOAT");
+    TCLAP::ValueArg<double> RES("", "cut_off", "Resolution cut-off (in A)", false, 1.0, "FLOAT");
+
+    INPUT_VOLUME.forceRequired();
+    INPUT_PAR.forceRequired();
+    OUTPUT_VOLUME.forceRequired();
+
+    //Add arguments
+    exe.add(OUTPUT_VOLUME);
+    exe.add(RES);
+    exe.add(NUM_PARTICLES);
+    exe.add(APIX);
+    exe.add(INPUT_PAR);
+    exe.add(INPUT_VOLUME);
+
+
+    //Parse the arguments
+    exe.parse(argc, argv);
+
+    //Prepare the input
+    Volume2DX input;
+    std::cout << "Reading the particles... \n";
+    input.read_volume(INPUT_VOLUME.getValue());
+
+    //Get number of particles to process
+    int number_particles = input.nz();
+    if(NUM_PARTICLES.isSet() && NUM_PARTICLES.getValue() < number_particles) {
+        number_particles = NUM_PARTICLES.getValue();
+    } 
+    std::cout << number_particles << " particles would be processed..\n";
+
+    //Get apix
+    double apix = APIX.getValue();
+    std::cout << "Pixel size is set to: " << apix << " A/px \n";
+
+    //Read the par file
+    std::vector<std::vector<double> > parameters(number_particles, std::vector<double>(5, 0.0));
+    tdx::File inFile (INPUT_PAR.getValue(), tdx::File::in);
+    if (!inFile.exists()){
+        std::cerr << "File not found: " << INPUT_PAR.getValue() << std::endl;
+        exit(1);
+    }
+    while(!inFile.eof()) {
+        std::string line = inFile.read_line();
+        if(line.length()==0) continue;
+        if(line[0] == 'C' || line[0] == 'c') continue;
+        std::vector<std::string> values = tdx::String::split(tdx::String::trim(line), ' ');
+        if(values.size() < 6) {
+            //std::cerr << "Omitting " << line << "\n";
+            continue;
+        }
+        int particle = std::stoi(values[0])-1;
+        if(particle < number_particles) {
+            parameters[particle] = {std::stod(values[1]), std::stod(values[2]), std::stod(values[3]), std::stod(values[4]), std::stod(values[5])};
+        }
+    }
+
+    tdx::data::MillerToPeakMultiMap all_reflections;
+    for(int i=0; i< number_particles; ++i) {
+        std::cout << "Processing particle: " << i+1 <<std::endl;
+        Volume2DX slice = input.get_slice(i);
+        double psi = parameters[i][0] * M_PI / 180;
+        double theta = parameters[i][1] * M_PI / 180;
+        double phi = parameters[i][2] * M_PI / 180;
+        double x_shift = parameters[i][3] / apix;
+        double y_shift = parameters[i][4] / apix;
+
+        std::cout << psi << " " << theta << " " << phi << " " << x_shift << " " << y_shift << "\n";
+
+        //Shift
+        slice.shift_volume(x_shift, y_shift, 0.0);
+
+        //Transform
+        auto data = slice.get_fourier();
+        for(const auto& itr : data) {
+            add_index(itr.first, itr.second.value(), all_reflections, psi, theta, phi);
+        }
+
+    }
+
+    std::cout << "Total number of reflections found in 3D Fourier space: " << all_reflections.size() << std::endl;
+    tdx::data::MillerToPeakMap averaged_data;
+    tdx::utilities::fourier_utilities::average_peaks(all_reflections, averaged_data);
+    tdx::data::ReflectionData fourier_data;
+    fourier_data.reset(averaged_data);
+
+    Volume2DX output(input.nx(), input.ny(), std::max(input.nx(), input.ny()));
+    output.set_fourier(fourier_data);
+    output.low_pass(RES.getValue()/apix);
+
+    output.write_volume(OUTPUT_VOLUME.getValue(), "mrc");
+
+    return 0;
+}

kernel/executables/2dx_processor.cpp

@@ -29,6 +29,9 @@ int main(int argc, char* argv[])
     exe.add(args::templates::BFACTOR);
     exe.add(args::templates::PSF);
     exe.add(args::templates::ZERO_PHASES);
+    exe.add(args::templates::SHIFTZ);
+    exe.add(args::templates::SHIFTY);
+    exe.add(args::templates::SHIFTX);
     exe.add(args::templates::INVERTZ);
     exe.add(args::templates::INVERTY);
     exe.add(args::templates::INVERTX);
@@ -119,6 +122,19 @@ int main(int argc, char* argv[])
     if(args::templates::INVERTY.getValue()) input.invert_hand(2);
     if(args::templates::INVERTZ.getValue()) input.invert_hand(3);
 
+
+    if(args::templates::SHIFTX.isSet() || args::templates::SHIFTY.isSet() ||  args::templates::SHIFTZ.isSet()) 
+    {
+        double x_shift = args::templates::SHIFTX.getValue();
+        double y_shift = args::templates::SHIFTY.getValue();
+        double z_shift = args::templates::SHIFTZ.getValue();
+
+        std::cout << "Shifting volume by: " << x_shift << " " << y_shift << " " << z_shift << "\n";
+
+        input.shift_volume(x_shift, y_shift, z_shift);
+
+    }
+
     if(args::templates::ZERO_PHASES.getValue())
     {
         std::cout << "Setting all phases to zero..\n";

kernel/toolkit/arguments/argument_template.hpp

@@ -48,6 +48,9 @@ namespace tdx
             static TCLAP::ValueArg<double> SLAB("", "slab", "The membrane height in ratio of the Z length of the volume", false, 1.0,"FLOAT");
             static TCLAP::ValueArg<std::string> TEMP_LOC("", "temp", "Folder to keep temporary files for each iteration (if not specified temp files will not be written)", false, "","FOLDER");
             static TCLAP::ValueArg<double> MASK_RES("", "mask-res", "The volume in each iteration will be lowpassed with this resolution, before generating a thresholded mask.", false, 15.0,"FLOAT");
+            static TCLAP::ValueArg<double> SHIFTX("", "x_shift", "The shift in x direction to be performed", false, 0.0,"FLOAT");
+            static TCLAP::ValueArg<double> SHIFTY("", "y_shift", "The shift in y direction to be performed", false, 0.0,"FLOAT");
+            static TCLAP::ValueArg<double> SHIFTZ("", "z_shift", "The shift in z direction to be performed", false, 0.0,"FLOAT");
 
             static TCLAP::SwitchArg INVERTED("", "inverted", "Produce an output map with inverted hand in all x,y,z direction", false);
             static TCLAP::SwitchArg INVERTX("", "invertx", "Produce an output map with inverted hand in x direction", false);

kernel/toolkit/data_structures/miller_index.cpp

@@ -50,11 +50,6 @@ bool ds::MillerIndex::operator <(const MillerIndex& rhs) const
           );
 }
 
-std::ostream& ds::MillerIndex::operator <<(std::ostream& os) const
-{
-    return os << h() << " " << k() << " " << l();
-}
-
 void ds::MillerIndex::set_values(int h, int k, int l)
 {
     this->initialize(h, k, l);

kernel/toolkit/data_structures/miller_index.hpp

@@ -71,7 +71,10 @@ namespace tdx
              * @param os
              * @return 
              */
-            std::ostream& operator<<(std::ostream& os) const;
+            friend inline std::ostream& operator<<(std::ostream& os, const MillerIndex& obj) {
+                os << " ( " << obj.h() << " " << obj.k() << " " << obj.l() << ") ";
+                return os;
+            }
 
             /**
              * Getter method for index h (x-dimension).

kernel/toolkit/data_structures/volume2dx.cpp

@@ -174,7 +174,7 @@ void ds::Volume2DX::set_real(const RealSpaceData& real)
 
 void ds::Volume2DX::fourier_from_real()
 {
-    std::cout << "Converting the real data to Fourier. \n";
+    //std::cout << "Converting the real data to Fourier. \n";
     if(_type == REAL)
     {
         _fourier.clear();
@@ -194,7 +194,7 @@ void ds::Volume2DX::fourier_from_real()
 
 void ds::Volume2DX::real_from_fourier()
 {
-    std::cout << "Converting the Fourier data to real. \n";
+    //std::cout << "Converting the Fourier data to real. \n";
     if(_type == FOURIER){
         double* real_data = fftw_alloc_real(nx()*ny()*nz());
         fftw_complex* complex_data = _fourier.fftw_data(fx(), fy(), fz());
@@ -569,6 +569,26 @@ ds::BinnedData ds::Volume2DX::calculate_structure_factors(double min_freq, doubl
     return binned_data;
 }
 
+void ds::Volume2DX::shift_volume(double x_shift, double y_shift, double z_shift) {
+    int dim_x = nx();
+    int dim_y = ny();
+    int dim_z = nz();
+
+    tdx::data::ReflectionData current, shifted;
+
+    current = get_fourier();
+    for (const auto& data : current) {
+        auto index = data.first;
+        double amp = data.second.amplitude();
+        double phase = data.second.phase()  - (2 * M_PI) * ((x_shift * index.h() * 1.0/dim_x)  + (y_shift * index.k() * 1.0/dim_y) + (z_shift * index.l() * 1.0/dim_z));
+        tdx::Complex shift_complex(amp*cos(phase), amp*sin(phase));
+        shifted.set_spot_at(index.h(), index.k(), index.l(), shift_complex, data.second.weight());
+    }
+
+    set_fourier(shifted);
+}
+
+
 ds::BinnedData ds::Volume2DX::fourier_shell_correlation(ds::Volume2DX reference, double min_freq, double max_freq, int resolution_bins)
 {
     BinnedData binnedFSC(min_freq, max_freq, resolution_bins);
@@ -1205,6 +1225,33 @@ ds::Volume2DX ds::Volume2DX::average2D(char axis)
     return averaged;
 }
 
+ds::Volume2DX ds::Volume2DX::get_slice(int slice_no) {
+    VolumeHeader head = header();
+    RealSpaceData current_data = get_real();
+    RealSpaceData new_data;
+
+    if(slice_no >= nz() || slice_no < 0) {
+        std::cerr << "Can't get slice: " << slice_no << " Range(0, " << nz() << ")\n";
+        exit(1);
+    }
+
+    head.set_mz(1);
+    head.set_sections(1);
+    new_data = RealSpaceData(nx(), ny(), 1);
+    for(int ix=0; ix<current_data.nx(); ix++)
+    {
+        for(int iy=0; iy<current_data.ny(); iy++)
+        {
+            new_data.set_value_at(ix, iy, 0, current_data.get_value_at(ix, iy, slice_no));
+        }
+    }
+
+    Volume2DX slice(head);
+    slice.set_real(new_data);
+    return slice;
+}
+
+
 ds::Volume2DX ds::Volume2DX::subsample(int factor)
 {
     std::cout << "Sub-sampling volume by " << factor << " times.\n";

kernel/toolkit/data_structures/volume2dx.hpp

@@ -350,6 +350,13 @@ namespace tdx
              */
             Volume2DX average2D(char axis);
 
+            /**
+             * Get the slice from the real space
+             * @param slice_no
+             * @return 
+             */
+            Volume2DX get_slice(int slice_no);
+
             /**
              * Fetches the real valued density at (x,y,z) from the volume
              * @param x
@@ -504,6 +511,12 @@ namespace tdx
              */
             Volume2DX generate_bead_model(int no_of_beads, double density_threshold, double bead_model_resolution = 2.0);
 
+            /*
+             * Shifts the volume to given values. Internally first gets fourier
+             * transform and then changes the phases.
+             */
+            void shift_volume(double x_shift, double y_shift, double z_shift);
+
             /**
              * Centers the density along the z axis. Internally adds PI*miller_index_l 
              * to the phase

kernel/toolkit/file_io/mrc_io.cpp

@@ -68,6 +68,10 @@ tdx::data::VolumeHeader io::mrc::get_header(const std::string file_name, const s
         header.set_ylen(infile.read_float());
         header.set_zlen(infile.read_float());
 
+        if(header.xlen() < 1.0) header.set_xlen(1.0);
+        if(header.ylen() < 1.0) header.set_ylen(1.0);
+        if(header.zlen() < 1.0) header.set_zlen(1.0);
+
         float alpha = infile.read_float();
         float beta = infile.read_float();
 

kernel/toolkit/transforms/fourier_transform_fftw.cpp

@@ -90,14 +90,14 @@ double ft::FourierTransformFFTW::NormalizationFactor()
 
 void ft::FourierTransformFFTW::Replan(double* real_data, fftw_complex* complex_data, int nx, int ny, int nz)
 {
-    std::cout << "Re-planning.. ";
+    //std::cout << "Re-planning.. ";
     _nx = nx;
     _ny = ny;
     _nz = nz;
     _plan_r2c = new fftw_plan(fftw_plan_dft_r2c_3d(nz, ny, nx, real_data, complex_data, FFTW_ESTIMATE));
     _plan_c2r = new fftw_plan(fftw_plan_dft_c2r_3d(nz, ny, nx, complex_data, real_data, FFTW_ESTIMATE));
     _plans_initialized = true;
-    std::cout << "Created new plans\n";
+    //std::cout << "Created new plans\n";
 }
 
 void ft::FourierTransformFFTW::RealToComplex(int nx, int ny, int nz, double* real_data, fftw_complex* complex_data)