Index: src/emlib/MRCReaderWriter.cpp =================================================================== --- src/emlib/MRCReaderWriter.cpp (revision 125) +++ src/emlib/MRCReaderWriter.cpp (working copy) @@ -1,4 +1,5 @@ #include "MRCReaderWriter.h" +#include // Reads an MRC file and translates the header to the general DensityHeader @@ -6,7 +7,7 @@ { // Read file filename.assign(fn_in); - read(data); + if (read(data)) return 1; // Translate header to DensityHeader header.ToDensityHeader(head); return 0; @@ -34,17 +35,20 @@ if (read_header()==1) { std::cout << "MRCReaderWriter::read >> Error reading MRC header of file " << filename << std::endl; - return 1; + throw std::ios_base::failure(""); } // Allocate memory size_t n = header.nx*header.ny*header.nz; + std::cout << "Found " << n << " voxels" << std::endl; (*pt)= new float[n]; // read if (read_data(*pt) == 1) { std::cout << "MRCReaderWriter::read >> Error reading MRC data of file " << filename << std::endl; - return 1; + throw std::ios_base::failure(""); } + } else { + throw std::ios_base::failure(""); } fs.close(); return 0; @@ -98,17 +102,19 @@ { int needswap,ierr; if (seek_to_data() != 0) - return 1; + throw std::ios_base::failure(""); size_t n=header.nx*header.ny*header.nz; // size of the grid ierr=read_grid(pt,sizeof(float), n); - if(ierr==1) - return ierr; // There is an error + if(ierr==1) { + throw std::ios_base::failure(""); // There is an error + } // Check for the necessity of changing the endian needswap = 0; for(size_t i=0;i 1e10) { + if (pt[i] > 1e10 || pt[i] < -1e10) { // Really big values usually result if the endian is not correct + std::cout << "Needing to swap order" << std::endl; needswap = 1; break; } @@ -186,13 +192,13 @@ if(!s.fail()) { // Write header if(write_header(s)==1) { - std::cout << "MRCReaderWriter::write. Error writing MRC header to file " + std::cerr << "MRCReaderWriter::write. Error writing MRC header to file " << fn << std::endl; return 1; } // Write values if(write_data(s,pt)==1) { - std::cout << "MRCReaderWriter::write. Error writing MRC data to file " + std::cerr << "MRCReaderWriter::write. Error writing MRC data to file " << fn << std::endl; return 1; } @@ -242,7 +248,7 @@ s.write((char *) &header.nlabl,wordsize); s.write((char *) &header.labels,sizeof(char)*MRC_NUM_LABELS*MRC_LABEL_SIZE); if(s.bad()) { - std::cout << "MRCReaderWriter::write_header. Error writing MRC header " + std::cerr << "MRCReaderWriter::write_header. Error writing MRC header " << "to file" << std::endl; return 1; } @@ -257,7 +263,7 @@ s.write((char *)pt,sizeof(float)*header.nx * header.ny * header.nz); if(s.bad()){ - std::cout << "MRCReaderWriter::write_header. Error writing MRC data to file" + std::cerr << "MRCReaderWriter::write_header. Error writing MRC data to file" << std::endl; return 1; } Index: src/emlib/DensityMap.cpp =================================================================== --- src/emlib/DensityMap.cpp (revision 125) +++ src/emlib/DensityMap.cpp (working copy) @@ -67,6 +67,9 @@ header_.nx = nx; header_.ny = ny; header_.nz = nz; + for (int i=0; i< nvox; ++i) { + data_[i]=0; + } } @@ -75,8 +78,8 @@ // TODO: we need to decide who does the allocation ( mapreaderwriter or // density)? if we keep the current implementation ( mapreaderwriter ) // we need to pass a pointer to data_ - std::cout<<"start" << std::endl; - float *f_data; + //std::cout<<"start" << std::endl; + float *f_data=0; if (reader.Read(filename, &f_data, header_) != 0) { std::cerr << " DensityMap::Read unable to read map encoded in file : " << filename << std::endl; @@ -87,9 +90,9 @@ normalized_ = false; calcRMS(); calc_all_voxel2loc(); - std::cout<<"before computer top" << std::endl; + //std::cout<<"before computer top" << std::endl; header_.compute_xyz_top(); - std::cout<<"after computer top" << std::endl; + //std::cout<<"after computer top" << std::endl; } void DensityMap::float2real(float *f_data, emreal **r_data) @@ -160,11 +163,59 @@ } } + +emreal DensityMap::get_value(unsigned int i, unsigned int j, unsigned int k) const { + if (i >= static_cast(header_.nx) + || j >= static_cast(header_.ny) + || k >= static_cast(header_.nz)) { + return 0; + } + return data_[k * header_.nx * header_.ny + j * header_.nx + i]; +} + emreal DensityMap::get_value(float x, float y, float z) const { return data_[loc2voxel(x,y,z)]; } +emreal DensityMap::get_smoothed_value(float x, float y, float z) const +{ + // trilirp in z, y, x + float c[3]={x,y,z}; + for (unsigned int i=0; i< 3; ++i){ + if (c[i] < header_.get_origin(i) || c[i] >= header_.get_top(i)) { + return 0; + } + } + + emreal fvox[3]; + unsigned int ivox[3]; + emreal r[3]; + for (unsigned int i=0; i< 3; ++i) { + fvox[i]= (c[i]- header_.get_origin(i))/header_.Objectpixelsize; + ivox[i]= static_cast(floor(fvox[i])); + r[i]= fvox[i]-ivox[i]; + assert(r[i] < 1.01 && r[i] >= -.01); + } + + emreal is[4]; + for (unsigned int i=0; i< 4; ++i) { + // operator >> has high precidence compared. Go fig. + unsigned int bx= ((i&2) >> 1); + unsigned int by= (i&1); + assert((bx==0 || bx==1) && (by==0 || by==1)); + is[i]= get_value(ivox[0]+bx, ivox[1]+by, ivox[2])*(1-r[2]) + + get_value(ivox[0]+bx, ivox[1]+by, ivox[2]+1U)*(r[2]); + } + + emreal js[2]; + for (unsigned int i=0; i< 2; ++i) { + js[i]= is[i*2] * (1- r[1]) + is[i*2+1] * (r[1]); + } + + return js[0]*(1-r[0]) + js[1]*(r[0]); +} + void DensityMap::calc_all_voxel2loc() { if (loc_calculated_) Index: src/emlib/DensityMap.h =================================================================== --- src/emlib/DensityMap.h (revision 125) +++ src/emlib/DensityMap.h (working copy) @@ -5,10 +5,12 @@ #include "EM_config.h" #include "DensityHeader.h" #include "Vector3.h" +#include "MapReaderWriter.h" + +#include #include #include #include -#include "MapReaderWriter.h" //! Class for handling density maps. /** CONVENTIONS @@ -83,6 +85,16 @@ */ emreal get_value(float x,float y,float z) const; + //! Gets the value at the coordinates interpolated from the center values + /** Trilinear interpolation is used. Any place outside the volume + has the value of 0. + */ + emreal get_smoothed_value(float x, float y, float z) const; + + //! Get the value of the ith, jth, kth voxel + /** return 0 if it is out of range */ + emreal get_value(unsigned int i, unsigned int j, unsigned int k) const; + //! Sets the origin of the header void set_origin(float x,float y,float z); Index: src/emlib/DensityHeader.h =================================================================== --- src/emlib/DensityHeader.h (revision 125) +++ src/emlib/DensityHeader.h (working copy) @@ -2,10 +2,12 @@ #define _DENSITYHEADER_H #include "EM_config.h" +#include "def.h" + #include #include #include -#include "def.h" +#include /** \todo change so that the att will not be encoded but loaded from @@ -132,6 +134,7 @@ float Exposuretime; //Exposuretime float Objectpixelsize; //this is the actual pixelsize float Microscope; //Microscope + //! Use this to set the size of the voxel float Pixelsize; //Pixelsize - used for the microscope CCD camera float CCDArea; //CCDArea float Defocus; //Defocus @@ -153,6 +156,14 @@ inline float get_xorigin() const {return xorigin;} inline float get_yorigin() const {return yorigin;} inline float get_zorigin() const {return zorigin;} + inline float get_origin(int i) const { + assert(i>=0 && i < 3); + switch (i) { + case 0: return xorigin; + case 1: return yorigin; + default: return zorigin; + } + } inline void set_xorigin(float x) {xorigin=x; top_calculated=false;} inline void set_yorigin(float y) {yorigin=y; top_calculated=false;} inline void set_zorigin(float z) {zorigin=z; top_calculated=false;}