00001
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036 #include <cstring>
00037 #include "fitsio.h"
00038 #include "portable_fileio.h"
00039 #include "geometry.h"
00040 #include "util.h"
00041 #include "ctf.h"
00042
00043 using namespace EMAN;
00044
00045 FitsIO::FitsIO(const string & fits_filename, IOMode rw)
00046 : filename(fits_filename), rw_mode(rw)
00047 {
00048 is_big_endian = ByteOrder::is_host_big_endian();
00049 is_new_file = false;
00050 initialized = false;
00051 fitsfile=0;
00052 }
00053
00054 FitsIO::~FitsIO()
00055 {
00056 if (fitsfile) {
00057 fclose(fitsfile);
00058 fitsfile = 0;
00059 }
00060 }
00061
00062 void FitsIO::init()
00063 {
00064 ENTERFUNC;
00065
00066 if (initialized) {
00067 return;
00068 }
00069
00070 initialized = true;
00071 fitsfile = sfopen(filename, rw_mode, &is_new_file);
00072
00073 EXITFUNC;
00074 }
00075
00076
00077 bool FitsIO::is_image_big_endian()
00078 {
00079 init();
00080 return is_big_endian;
00081 }
00082
00083 bool FitsIO::is_valid(const void *first_block, off_t)
00084 {
00085 ENTERFUNC;
00086
00087 if (!first_block) {
00088 return false;
00089 }
00090
00091 if (strncmp("SIMPLE ",(const char *)first_block,8)==0) return true;
00092
00093 EXITFUNC;
00094 return false;
00095 }
00096
00097 int FitsIO::read_header(Dict & dict, int image_index, const Region * area, bool )
00098 {
00099 ENTERFUNC;
00100
00101
00102
00103
00104 if(image_index == -1) {
00105 image_index = 0;
00106 }
00107
00108 if(image_index != 0) {
00109 throw ImageReadException(filename, "no stack allowed for MRC image. For take 2D slice out of 3D image, read the 3D image first, then use get_clip().");
00110 }
00111 init();
00112
00113 if (area) throw ImageReadException(filename,"Area reading not supported for FITS format");
00114
00115 dict["nx"]=1;
00116 dict["ny"]=1;
00117 dict["nz"]=1;
00118
00119 char s[81],lbl[9],val[80];
00120 int dim=0;
00121 s[80]=0;
00122 rewind(fitsfile);
00123 for (fread(s,80,1,fitsfile); strncmp("END",s,3); fread(s,80,1,fitsfile)) {
00124 sscanf(s,"%8s = %[^/]",lbl,val);
00125
00126 if (strncmp("SIMPLE ",s,8)==0) continue;
00127 else if (strncmp("END ",s,8)==0) break;
00128
00129 else if (strncmp("NAXIS ",s,8)==0) dim=atoi(val);
00130 else if (strncmp("NAXIS",s,5)==0) {
00131 if (s[5]=='1') dict["nx"]=atoi(val);
00132 if (s[5]=='2') dict["ny"]=atoi(val);
00133 if (s[5]=='3') dict["nz"]=atoi(val);
00134 }
00135 else {
00136 dict[(string)"FITS."+lbl]=val;
00137 }
00138 }
00139
00140 dstart=((ftell(fitsfile)-1)/2880+1)*2880;
00141
00142 int xlen = 0, ylen = 0, zlen = 0;
00143 dtype=atoi(dict["FITS.BITPIX"]);
00144 EMUtil::get_region_dims(area, dict["nx"], &xlen, dict["ny"], &ylen, dict["nz"], &zlen);
00145
00146 dict["nx"] = nx=xlen;
00147 dict["ny"] = ny=ylen;
00148 dict["nz"] = nz=zlen;
00149
00150 EXITFUNC;
00151 return 0;
00152 }
00153
00154 int FitsIO::write_header(const Dict &, int, const Region*, EMUtil::EMDataType, bool)
00155 {
00156 ENTERFUNC;
00157
00158 LOGWARN("FITS write is not supported.");
00159 EXITFUNC;
00160 return 0;
00161 }
00162
00163 int FitsIO::read_data(float *rdata, int image_index, const Region *, bool )
00164 {
00165 ENTERFUNC;
00166 size_t i;
00167 size_t size = (size_t)nx*ny*nz;
00168
00169
00170 image_index = 0;
00171 check_read_access(image_index, rdata);
00172
00173 portable_fseek(fitsfile, dstart, SEEK_SET);
00174 char *cdata=(char *)rdata;
00175 short *sdata=(short *)rdata;
00176 int *idata=(int *)rdata;
00177 double *ddata;
00178
00179 switch (dtype) {
00180 case 8:
00181 fread(cdata,nx,ny*nz,fitsfile);
00182 for (i=size-1; i<size; i--) rdata[i]=cdata[i];
00183 break;
00184 case 16:
00185 fread(cdata,nx,ny*nz*2,fitsfile);
00186 if (!ByteOrder::is_host_big_endian()) ByteOrder::swap_bytes((short*) sdata, size);
00187 for (i=size-1; i<size; i--) rdata[i]=sdata[i];
00188 break;
00189 case 32:
00190 fread(cdata,nx,ny*nz*4,fitsfile);
00191 if (!ByteOrder::is_host_big_endian()) ByteOrder::swap_bytes((int*) rdata, size);
00192 for (i=0; i<size; i++) rdata[i]=static_cast<float>(idata[i]);
00193 break;
00194 case -32:
00195 fread(cdata,nx*4,ny*nz,fitsfile);
00196 if (!ByteOrder::is_host_big_endian()) ByteOrder::swap_bytes((float*) rdata, size);
00197 break;
00198 case -64:
00199 ddata=(double *)malloc(size*8);
00200 fread(ddata,nx,ny*nz*8,fitsfile);
00201 if (!ByteOrder::is_host_big_endian()) ByteOrder::swap_bytes((double*) ddata, size);
00202 for (i=0; i<size; i++) rdata[i]=static_cast<float>(ddata[i]);
00203 free(ddata);
00204 break;
00205 }
00206
00207 EXITFUNC;
00208 return 0;
00209 }
00210
00211 int FitsIO::write_data(float *data, int image_index, const Region*,
00212 EMUtil::EMDataType, bool)
00213 {
00214 ENTERFUNC;
00215
00216 check_write_access(rw_mode, image_index, 1, data);
00217
00218
00219
00220 EXITFUNC;
00221 return 0;
00222 }
00223
00224
00225 bool FitsIO::is_complex_mode()
00226 {
00227 init();
00228 return false;
00229 }
00230
00231 void FitsIO::flush()
00232 {
00233 fflush(fitsfile);
00234 }
00235
00236 int FitsIO::read_ctf(Ctf &, int)
00237 {
00238 ENTERFUNC;
00239 init();
00240 EXITFUNC;
00241 return 0;
00242 }
00243
00244 void FitsIO::write_ctf(const Ctf &, int)
00245 {
00246 ENTERFUNC;
00247 init();
00248
00249 EXITFUNC;
00250 }