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 <exception>
// 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<n;i++)
- if (pt[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<unsigned int>(header_.nx)
+ || j >= static_cast<unsigned int>(header_.ny)
+ || k >= static_cast<unsigned int>(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<int>(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 <cassert>
#include <iostream>
#include <iomanip>
#include <stdexcept>
-#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 <iostream>
#include <fstream>
#include <cstring>
-#include "def.h"
+#include <cassert>
/** \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;}