EMAN2
omapio.cpp
Go to the documentation of this file.
00001 
00005 /*
00006  * Author: Grant Tang, 06/07/2011 (gtang@bcm.edu)
00007  * Copyright (c) 2000-2006 Baylor College of Medicine
00008  *
00009  * This software is issued under a joint BSD/GNU license. You may use the
00010  * source code in this file under either license. However, note that the
00011  * complete EMAN2 and SPARX software packages have some GPL dependencies,
00012  * so you are responsible for compliance with the licenses of these packages
00013  * if you opt to use BSD licensing. The warranty disclaimer below holds
00014  * in either instance.
00015  *
00016  * This complete copyright notice must be included in any revised version of the
00017  * source code. Additional authorship citations may be added, but existing
00018  * author citations must be preserved.
00019  *
00020  * This program is free software; you can redistribute it and/or modify
00021  * it under the terms of the GNU General Public License as published by
00022  * the Free Software Foundation; either version 2 of the License, or
00023  * (at your option) any later version.
00024  *
00025  * This program is distributed in the hope that it will be useful,
00026  * but WITHOUT ANY WARRANTY; without even the implied warranty of
00027  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
00028  * GNU General Public License for more details.
00029  *
00030  * You should have received a copy of the GNU General Public License
00031  * along with this program; if not, write to the Free Software
00032  * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
00033  *
00034  * */
00035 
00036 #include "omapio.h"
00037 #include "portable_fileio.h"
00038 
00039 using namespace EMAN;
00040 
00041 OmapIO::OmapIO(const string & omapname, IOMode rw) :
00042                 filename(omapname), rw_mode(rw), omapfile(0),
00043                 is_big_endian(false), initialized(false), is_new_file(false)
00044 {
00045         memset(&omaph, 0, sizeof(OmapHeader));
00046         is_big_endian = ByteOrder::is_host_big_endian();
00047 }
00048 
00049 OmapIO::~OmapIO()
00050 {
00051         if (omapfile) {
00052                 fclose(omapfile);
00053                 omapfile = 0;
00054         }
00055 }
00056 
00057 void OmapIO::init()
00058 {
00059         ENTERFUNC;
00060 
00061         if (initialized) {
00062                 return;
00063         }
00064 
00065         initialized = true;
00066         omapfile = sfopen(filename, rw_mode, &is_new_file);
00067 
00068         char record[512];       //record size 512 bytes
00069 
00070         if (!is_new_file) {
00071                 if (fread(record, 512, 1, omapfile) != 1) {
00072                         throw ImageReadException(filename, "OMAP header");
00073                 }
00074 
00075                 for(int i=0; i<512; i++) {
00076                         if(!isprint(record[i])) {
00077                                 portable_fseek(omapfile, 0, SEEK_SET);
00078                                 break;
00079                         }
00080 
00081                         if(record[i] == '\0') break;
00082                 }
00083 
00084                 if (fread(&omaph, sizeof(OmapHeader), 1, omapfile) != 1) {
00085                         throw ImageReadException(filename, "OMAP header");
00086                 }
00087 
00088                 if (!is_valid(&omaph)) {
00089                         throw ImageReadException(filename, "invalid OMAP");
00090                 }
00091 
00092                 if(!ByteOrder::is_host_big_endian()) {  //omap first record is always  big endian
00093                         ByteOrder::swap_bytes((short *) &omaph, 256);   //each record 512 bytes, 256 short intergers
00094                 }
00095         }
00096 
00097         EXITFUNC;
00098 }
00099 
00100 int OmapIO::read_header(EMAN::Dict& dict, int, EMAN::Region const*, bool)
00101 {
00102         ENTERFUNC;
00103         init();
00104 
00105         dict["OMAP.xstart"] = (int)omaph.xstart;
00106         dict["OMAP.ystart"] = (int)omaph.ystart;
00107         dict["OMAP.zstart"] = (int)omaph.zstart;
00108         dict["nx"] = (int)omaph.nx;
00109         dict["ny"] = (int)omaph.ny;
00110         dict["nz"] = (int)omaph.nz;
00111         dict["apix_x"] = (int)omaph.apix_x;
00112         dict["apix_y"] = (int)omaph.apix_y;
00113         dict["apix_z"] = (int)omaph.apix_z;
00114         dict["OMAP.cellA"] = (int)(omaph.header10/omaph.scale);
00115         dict["OMAP.cellB"] = (int)(omaph.header11/omaph.scale);
00116         dict["OMAP.cellC"] = (int)(omaph.header12/omaph.scale);
00117         dict["alpha"] = (int)(omaph.alpha/omaph.scale);
00118         dict["beta"] = (int)(omaph.beta/omaph.scale);
00119         dict["gamma"] = (int)(omaph.gamma/omaph.scale);
00120         dict["OMAP.iprod"] = (int)omaph.iprod;
00121         dict["OMAP.iplus"] = (int)omaph.iplus;
00122         dict["OMAP.scale"] = (int)omaph.scale?omaph.scale:100;
00123         dict["OMAP.scale2"] = (int)omaph.scale2?omaph.scale2:100;
00124 
00125         float prod = (float)(omaph.iprod/(int)dict["OMAP.scale2"]);
00126         float plus = (float)omaph.iplus;
00127 
00128         dict["OMAP.min"] = (omaph.imin - plus)/prod;
00129         dict["OMAP.imax"] = (omaph.imax - plus)/prod;
00130         dict["OMAP.sigma"] = (omaph.isigma - plus)/prod;
00131         dict["OMAP.mean"] = (omaph.imean - plus)/prod;
00132 
00133         if((float)dict["OMAP.sigma"] < 0.001f || (float)dict["OMAP.sigma"] > 10.0f) {
00134                 std::cout << "Warning : Suspect value of sigma : " << (float)dict["OMAP.sigma"] << std::endl;
00135                 dict["OMAP.sigma"] = 0;         //flag bad sigma
00136         }
00137 
00138         EXITFUNC;
00139         return 0;
00140 }
00141 
00142 int OmapIO::read_data(float *rdata, int, EMAN::Region const*, bool)
00143 {
00144         ENTERFUNC;
00145 //      std::cout << "file pointer location = " << ftell(omapfile) << std::endl;
00146 
00147         // cubes in each direction
00148         int inx = omaph.nx/8;
00149         int iny = omaph.ny/8;
00150         int inz = omaph.nz/8;
00151 
00152         //include partial cube
00153         if(omaph.nx%8 > 0) ++inx;
00154         if(omaph.ny%8 > 0) ++iny;
00155         if(omaph.nz%8 > 0) ++inz;
00156 
00157         // How much element of last cube to read ?
00158         int xtraX = omaph.nx%8;
00159         int xtraY = omaph.ny%8;
00160         int xtraZ = omaph.nz%8;
00161 
00162 //      std::cout << "Total records = " << inx*iny*inz << std::endl;
00163 
00164         float prod = (float)omaph.iprod/omaph.scale2?omaph.scale2:100;
00165         float plus = omaph.iplus;
00166         float pixel = 0.0f;
00167         unsigned char record[512];
00168         for (int k=0; k < inz; k++) {
00169                 for (int j=0;  j < iny; j++) {
00170                         for (int i=0; i < inx; i++) {
00171                                 if (fread(record, 512, 1, omapfile) != 1) {
00172                                         throw ImageReadException(filename, "OMAP data");
00173                                 }
00174 
00175                                 //swap bytes for little endian
00176                                 bool byteswap = false;
00177                                 if(!ByteOrder::is_host_big_endian()) {
00178                                         byteswap = true;
00179                                 }
00180                                 if(byteswap) {
00181                                         for (int ii=0; ii < 511; ii+=2) {
00182                                                 char tempchar = record[ii];
00183                                                 record[ii] = record[ii+1];
00184                                                 record[ii+1] = tempchar;
00185                                         }
00186                                 }
00187 
00188                                 int cubieSizeX = 8;
00189                                 int cubieSizeY = 8;
00190                                 int cubieSizeZ = 8;
00191 
00192                                 //for partial cube
00193                                 if (xtraX > 0) if (i == inx-1) cubieSizeX = xtraX;
00194                                 if (xtraY > 0) if (j == iny-1) cubieSizeY = xtraY;
00195                                 if (xtraZ > 0) if (k == inz-1) cubieSizeZ = xtraZ;
00196 
00197                                 for (int n=0; n < cubieSizeZ; n++) {
00198                                         for (int m=0;  m < cubieSizeY; m++) {
00199                                                 for (int l=0; l < cubieSizeX; l++) {
00200                                                         unsigned char sboxLMN = record[8*8*n+8*m+l];
00201 
00202                                                         pixel = ((float)sboxLMN - plus)/prod;
00203 
00204                                                         int pt3 = k*8 + n;
00205                                                         int pt2 = j*8 + m;
00206                                                         int pt1 = i*8 + l;
00207 
00208                                                         rdata[pt3*omaph.nx*omaph.ny + pt2*omaph.nx + pt1] = (pixel-plus)/prod;
00209                                                         if(omaph.isigma>0) rdata[pt3*omaph.nx*omaph.ny + pt2*omaph.nx + pt1] /= (float)omaph.isigma;
00210                                                 }
00211                                         }
00212                                 }
00213 
00214                         }
00215                 }
00216         }
00217 
00218         EXITFUNC;
00219         return 0;
00220 }
00221 
00222 int OmapIO::write_header(EMAN::Dict const&, int, EMAN::Region const*, EMAN::EMUtil::EMDataType, bool)
00223 {
00224         ENTERFUNC;
00225 
00226         throw ImageWriteException("N/A", "No writing for Omap images");
00227 
00228         EXITFUNC;
00229         return 0;
00230 }
00231 
00232 int OmapIO::write_data(float*, int, EMAN::Region const*, EMAN::EMUtil::EMDataType, bool)
00233 {
00234         ENTERFUNC;
00235 
00236         EXITFUNC;
00237         return 0;
00238 }
00239 
00240 bool OmapIO::is_valid(const void *first_block, off_t file_size)
00241 {
00242         ENTERFUNC;
00243 
00244         if (!first_block) {
00245                 return false;
00246         }
00247 
00248         const short *data = static_cast < const short *>(first_block);
00249         short xstart = data[0];
00250         short ystart = data[1];
00251         short zstart = data[2];
00252         short nx = data[3];
00253         short ny = data[4];
00254         short nz = data[5];
00255         short const_value = data[18];
00256 
00257         if(!ByteOrder::is_host_big_endian()) {
00258                 ByteOrder::swap_bytes(&xstart);
00259                 ByteOrder::swap_bytes(&ystart);
00260                 ByteOrder::swap_bytes(&zstart);
00261                 ByteOrder::swap_bytes(&nx);
00262                 ByteOrder::swap_bytes(&ny);
00263                 ByteOrder::swap_bytes(&nz);
00264                 ByteOrder::swap_bytes(&const_value);
00265         }
00266 
00267         if(const_value != 100) return false;
00268         if(nx<=0 || ny<=0 || nz<=0 || nx>10000 || ny>10000 || nz>10000) return false;
00269 
00270         EXITFUNC;
00271         return true;
00272 }
00273 
00274 bool OmapIO::is_image_big_endian()
00275 {
00276         return true;
00277 }
00278 
00279 bool OmapIO::is_complex_mode()
00280 {
00281         return false;
00282 }
00283 
00284 void OmapIO::flush()
00285 {
00286         fflush(omapfile);
00287 }