ecat63w.c
Go to the documentation of this file.
00001 /******************************************************************************
00002 
00003   ecat63w.c  (c) 2003-2008 by Turku PET Centre
00004 
00005   Procedures for writing ECAT 6.3 matrix data.
00006 
00007   This program is free software; you can redistribute it and/or modify it under
00008   the terms of the GNU General Public License as published by the Free Software
00009   Foundation; either version 2 of the License, or (at your option) any later
00010   version.
00011 
00012   This program is distributed in the hope that it will be useful, but WITHOUT
00013   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00014   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
00015 
00016   You should have received a copy of the GNU General Public License along with
00017   this program; if not, write to the Free Software Foundation, Inc., 59 Temple
00018   Place, Suite 330, Boston, MA 02111-1307 USA.
00019 
00020   Turku PET Centre hereby disclaims all copyright interest in the program.
00021   Juhani Knuuti
00022   Director, Professor
00023   Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
00024 
00025   Assumptions:
00026   1. All data is always saved in little endian byte order (i386 and VAX).
00027   2. Data is automatically saved in one of the little endian formats
00028      as specified in header data_type.
00029   3. VAX data can be saved correctly only in 2-byte formats.
00030 
00031   Modification history:
00032   2003-07-21 Vesa Oikonen
00033     Created based on previous contents of ecat63.c.
00034   2003-09-08 VO
00035     Added function ecat63WriteImageMatrix() and ecat63WriteScanMatrix().
00036   2003-11-30 VO
00037     For now, calls temp_roundf() instead of roundf().
00038   2004-01-07 VO
00039     ecat63WriteImageMatrix(): corrected img min&max in header.
00040   2004-09-17 VO
00041     Doxygen style comments.
00042   2004-12-28 VO
00043     Included function ecat63_is_scaling_needed().
00044     This function is applied to determine whether scal factor can be set to
00045     one in case that all pixel values are close to integers and small enough.
00046   2007-01-27 VO
00047     Unsigned char pointer was corrected to signed in ecat63WriteMatdata().
00048   2007-09-02 VO
00049     Backup file extension changed from % to .bak.
00050 
00051 ******************************************************************************/
00052 #include <stdio.h>
00053 #include <stdlib.h>
00054 #include <math.h>
00055 #include <string.h>
00056 #include <unistd.h>
00057 #include <time.h>
00058 /*****************************************************************************/
00059 #include <swap.h>
00060 #include <petc99.h>
00061 #include "include/ecat63.h"
00062 /*****************************************************************************/
00063 
00064 /*****************************************************************************/
00073 int ecat63WriteMainheader(FILE *fp, ECAT63_mainheader *h) {
00074   char buf[MatBLKSIZE];
00075   int i, little, tovax;
00076 
00077 
00078   if(ECAT63_TEST) printf("ecat63WriteMainheader()\n");
00079   little=little_endian();
00080   /* Clear buf */
00081   memset(buf, 0, MatBLKSIZE);
00082   /* Check arguments */
00083   if(fp==NULL || h->data_type<1 || h->data_type>7) return(1);
00084   if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
00085     tovax=1; else tovax=0;
00086 
00087   /* Copy short ints to buf */
00088   memcpy(buf+50, &h->data_type, 2); memcpy(buf+48, &h->sw_version, 2);
00089   memcpy(buf+52, &h->system_type, 2); memcpy(buf+54, &h->file_type, 2);
00090   memcpy(buf+66, &h->scan_start_day, 2); memcpy(buf+68, &h->scan_start_month, 2);
00091   memcpy(buf+70, &h->scan_start_year, 2); memcpy(buf+72, &h->scan_start_hour, 2);
00092   memcpy(buf+74, &h->scan_start_minute, 2); memcpy(buf+76, &h->scan_start_second, 2);
00093   memcpy(buf+134, &h->rot_source_speed, 2); memcpy(buf+136, &h->wobble_speed, 2);
00094   memcpy(buf+138, &h->transm_source_type, 2); memcpy(buf+148, &h->transaxial_samp_mode, 2);
00095   memcpy(buf+150, &h->coin_samp_mode, 2); memcpy(buf+152, &h->axial_samp_mode, 2);
00096   memcpy(buf+158, &h->calibration_units, 2); memcpy(buf+160, &h->compression_code, 2);
00097   memcpy(buf+350, &h->acquisition_type, 2); memcpy(buf+352, &h->bed_type, 2);
00098   memcpy(buf+354, &h->septa_type, 2); memcpy(buf+376, &h->num_planes, 2);
00099   memcpy(buf+378, &h->num_frames, 2); memcpy(buf+380, &h->num_gates, 2);
00100   memcpy(buf+382, &h->num_bed_pos, 2); memcpy(buf+452, &h->lwr_sctr_thres, 2);
00101   memcpy(buf+454, &h->lwr_true_thres, 2); memcpy(buf+456, &h->upr_true_thres, 2);
00102   memcpy(buf+472, h->fill2, 40);
00103   /* big to little endian if necessary */
00104   if(!little) swabip(buf, MatBLKSIZE);
00105 
00106   /* Copy floats to buf */
00107   ecat63wFloat(&h->isotope_halflife, buf+86, tovax, little);
00108   ecat63wFloat(&h->gantry_tilt, buf+122, tovax, little);
00109   ecat63wFloat(&h->gantry_rotation, buf+126, tovax, little);
00110   ecat63wFloat(&h->bed_elevation, buf+130, tovax, little);
00111   ecat63wFloat(&h->axial_fov, buf+140, tovax, little);
00112   ecat63wFloat(&h->transaxial_fov, buf+144, tovax, little);
00113   ecat63wFloat(&h->calibration_factor, buf+154, tovax, little);
00114   ecat63wFloat(&h->init_bed_position, buf+384, tovax, little);
00115   for(i=0; i<15; i++) ecat63wFloat(&h->bed_offset[i], buf+388+4*i, tovax, little);
00116   ecat63wFloat(&h->plane_separation, buf+448, tovax, little);
00117   ecat63wFloat(&h->init_bed_position, buf+458, tovax, little);
00118 
00119   /* Copy chars */
00120   /*memcpy(buf+0, h->ecat_format, 14);*/
00121   memcpy(buf+14, h->fill1, 14);
00122   memcpy(buf+28, h->original_file_name, 20); memcpy(buf+56, h->node_id, 10);
00123   memcpy(buf+78, h->isotope_code, 8); memcpy(buf+90, h->radiopharmaceutical, 32);
00124   memcpy(buf+162, h->study_name, 12); memcpy(buf+174, h->patient_id, 16);
00125   memcpy(buf+190, h->patient_name, 32); buf[222]=h->patient_sex;
00126   memcpy(buf+223, h->patient_age, 10); memcpy(buf+233, h->patient_height, 10);
00127   memcpy(buf+243, h->patient_weight, 10); buf[253]=h->patient_dexterity;
00128   memcpy(buf+254, h->physician_name, 32); memcpy(buf+286, h->operator_name, 32);
00129   memcpy(buf+318, h->study_description, 32); memcpy(buf+356, h->facility_name, 20);
00130   memcpy(buf+462, h->user_process_code, 10);
00131 
00132   /* Write main header */
00133   fseek(fp, 0*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=0*MatBLKSIZE) return(2);
00134   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
00135 
00136   return(0);
00137 }
00138 /*****************************************************************************/
00139 
00140 /*****************************************************************************/
00150 int ecat63WriteImageheader(FILE *fp, int block, ECAT63_imageheader *h) {
00151   char buf[MatBLKSIZE];
00152   int i, little, tovax;
00153 
00154 
00155   if(ECAT63_TEST) printf("ecat63WriteImageheader(fp, %d, ih)\n", block);
00156   little=little_endian();
00157   /* Clear buf */
00158   memset(buf, 0, MatBLKSIZE);
00159   /* Check arguments */
00160   if(fp==NULL || block<3 || h->data_type<1 || h->data_type>7) return(1);
00161   if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
00162     tovax=1; else tovax=0;
00163 
00164   /* Copy short ints to buf */
00165   memcpy(buf+126, &h->data_type, 2); memcpy(buf+128, &h->num_dimensions, 2);
00166   memcpy(buf+132, &h->dimension_1, 2); memcpy(buf+134, &h->dimension_2, 2);
00167   memcpy(buf+176, &h->image_min, 2); memcpy(buf+178, &h->image_max, 2);
00168   memcpy(buf+200, &h->slice_location, 2); memcpy(buf+202, &h->recon_start_hour, 2);
00169   memcpy(buf+204, &h->recon_start_min, 2); memcpy(buf+206, &h->recon_start_sec, 2);
00170   memcpy(buf+236, &h->filter_code, 2); memcpy(buf+376, &h->processing_code, 2);
00171   memcpy(buf+380, &h->quant_units, 2); memcpy(buf+382, &h->recon_start_day, 2);
00172   memcpy(buf+384, &h->recon_start_month, 2); memcpy(buf+386, &h->recon_start_year, 2);
00173   memcpy(buf+460, h->fill2, 52);
00174   /* big to little endian if necessary */
00175   if(!little) swabip(buf, MatBLKSIZE);
00176 
00177   /* Copy floats to buf */
00178   ecat63wFloat(&h->x_origin, buf+160, tovax, little);
00179   ecat63wFloat(&h->y_origin, buf+164, tovax, little);
00180   ecat63wFloat(&h->recon_scale, buf+168, tovax, little);
00181   ecat63wFloat(&h->quant_scale, buf+172, tovax, little);
00182   ecat63wFloat(&h->pixel_size, buf+184, tovax, little);
00183   ecat63wFloat(&h->slice_width, buf+188, tovax, little);
00184   ecat63wFloat(&h->image_rotation, buf+296, tovax, little);
00185   ecat63wFloat(&h->plane_eff_corr_fctr, buf+300, tovax, little);
00186   ecat63wFloat(&h->decay_corr_fctr, buf+304, tovax, little);
00187   ecat63wFloat(&h->loss_corr_fctr, buf+308, tovax, little);
00188   ecat63wFloat(&h->intrinsic_tilt, buf+312, tovax, little);
00189   ecat63wFloat(&h->ecat_calibration_fctr, buf+388, tovax, little);
00190   ecat63wFloat(&h->well_counter_cal_fctr, buf+392, tovax, little);
00191   for(i=0; i<6; i++) ecat63wFloat(&h->filter_params[i], buf+396+4*i, tovax, little);
00192 
00193   /* Copy ints to buf */
00194   ecat63wInt(&h->frame_duration, buf+192, tovax, little);
00195   ecat63wInt(&h->frame_start_time, buf+196, tovax, little);
00196   ecat63wInt(&h->scan_matrix_num, buf+238, tovax, little);
00197   ecat63wInt(&h->norm_matrix_num, buf+242, tovax, little);
00198   ecat63wInt(&h->atten_cor_mat_num, buf+246, tovax, little);
00199 
00200   /* Copy chars */
00201   memcpy(buf+0, h->fill1, 126); memcpy(buf+420, h->annotation, 40);
00202 
00203   /* Write subheader */
00204   fseek(fp, (block-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(block-1)*MatBLKSIZE) return(2);
00205   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
00206 
00207   return(0);
00208 }
00209 /*****************************************************************************/
00210 
00211 /*****************************************************************************/
00221 int ecat63WriteAttnheader(FILE *fp, int block, ECAT63_attnheader *h) {
00222   unsigned char buf[MatBLKSIZE];
00223   int little, tovax;
00224 
00225   if(ECAT63_TEST) printf("ecat63WriteAttnheader(fp, %d, ah)\n", block);
00226   little=little_endian();
00227   /* Clear buf */
00228   memset(buf, 0, MatBLKSIZE);
00229   /* Check arguments */
00230   if(fp==NULL || block<3 || h->data_type<1 || h->data_type>7) return(1);
00231   if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
00232     tovax=1; else tovax=0;
00233 
00234   /* Copy short ints to buf */
00235   memcpy(buf+126, &h->data_type, 2); memcpy(buf+128, &h->attenuation_type, 2);
00236   memcpy(buf+132, &h->dimension_1, 2); memcpy(buf+134, &h->dimension_2, 2);
00237   /* big to little endian if necessary */
00238   if(!little) swabip(buf, MatBLKSIZE);
00239 
00240   /* Copy floats to buf */
00241   ecat63wFloat(&h->scale_factor, buf+182, tovax, little);
00242   ecat63wFloat(&h->x_origin, buf+186, tovax, little);
00243   ecat63wFloat(&h->y_origin, buf+190, tovax, little);
00244   ecat63wFloat(&h->x_radius, buf+194, tovax, little);
00245   ecat63wFloat(&h->y_radius, buf+198, tovax, little);
00246   ecat63wFloat(&h->tilt_angle, buf+202, tovax, little);
00247   ecat63wFloat(&h->attenuation_coeff, buf+206, tovax, little);
00248   ecat63wFloat(&h->sample_distance, buf+210, tovax, little);
00249 
00250   /* Write subheader */
00251   fseek(fp, (block-1)*MatBLKSIZE, SEEK_SET);
00252   if(ftell(fp)!=(block-1)*MatBLKSIZE) return(2);
00253   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
00254 
00255   return(0);
00256 }
00257 /*****************************************************************************/
00258 
00259 /*****************************************************************************/
00269 int ecat63WriteScanheader(FILE *fp, int block, ECAT63_scanheader *h) {
00270   unsigned char buf[MatBLKSIZE];
00271   int i, little, tovax;
00272 
00273 
00274   if(ECAT63_TEST) printf("ecat63WriteScanheader(fp, %d, ih)\n", block);
00275   little=little_endian();
00276   /* Clear buf */
00277   memset(buf, 0, MatBLKSIZE);
00278   /* Check arguments */
00279   if(fp==NULL || block<3 || h->data_type<1 || h->data_type>7) return(1);
00280   if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
00281     tovax=1; else tovax=0;
00282 
00283   /* Copy short ints to buf */
00284   memcpy(buf+126, &h->data_type, 2);
00285   memcpy(buf+132, &h->dimension_1, 2); memcpy(buf+134, &h->dimension_2, 2);
00286   memcpy(buf+136, &h->smoothing, 2); memcpy(buf+138, &h->processing_code, 2);
00287   memcpy(buf+170, &h->frame_duration_sec, 2);
00288   memcpy(buf+192, &h->scan_min, 2); memcpy(buf+194, &h->scan_max, 2);
00289   memcpy(buf+468, h->fill2, 44);
00290   /* big to little endian if necessary */
00291   if(!little) swabip(buf, MatBLKSIZE);
00292 
00293   /* Copy floats to buf */
00294   ecat63wFloat(&h->sample_distance, buf+146, tovax, little);
00295   ecat63wFloat(&h->isotope_halflife, buf+166, tovax, little);
00296   ecat63wFloat(&h->scale_factor, buf+182, tovax, little);
00297   for(i=0; i<16; i++) {
00298     ecat63wFloat(&h->cor_singles[i], buf+316+4*i, tovax, little);
00299     ecat63wFloat(&h->uncor_singles[i], buf+380+4*i, tovax, little);
00300   }
00301   ecat63wFloat(&h->tot_avg_cor, buf+444, tovax, little);
00302   ecat63wFloat(&h->tot_avg_uncor, buf+448, tovax, little);
00303   ecat63wFloat(&h->loss_correction_fctr, buf+464, tovax, little);
00304 
00305   /* Copy ints to buf */
00306   ecat63wInt(&h->gate_duration, buf+172, tovax, little);
00307   ecat63wInt(&h->r_wave_offset, buf+176, tovax, little);
00308   ecat63wInt(&h->prompts, buf+196, tovax, little);
00309   ecat63wInt(&h->delayed, buf+200, tovax, little);
00310   ecat63wInt(&h->multiples, buf+204, tovax, little);
00311   ecat63wInt(&h->net_trues, buf+208, tovax, little);
00312   ecat63wInt(&h->total_coin_rate, buf+452, tovax, little);
00313   ecat63wInt(&h->frame_start_time, buf+456, tovax, little);
00314   ecat63wInt(&h->frame_duration, buf+460, tovax, little);
00315 
00316   /* Copy chars */
00317   memcpy(buf+0, h->fill1, 126);
00318 
00319   /* Write subheader */
00320   fseek(fp, (block-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(block-1)*MatBLKSIZE) return(2);
00321   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
00322 
00323   return(0);
00324 }/*****************************************************************************/
00325 
00326 /*****************************************************************************/
00336 int ecat63WriteNormheader(FILE *fp, int block, ECAT63_normheader *h)
00337 {
00338   unsigned char buf[MatBLKSIZE];
00339   int little, tovax;
00340 
00341   if(ECAT63_TEST) printf("ecat63WriteNormheader(fp, %d, nh)\n", block);
00342   little=little_endian();
00343   /* Clear buf */
00344   memset(buf, 0, MatBLKSIZE);
00345   /* Check arguments */
00346   if(fp==NULL || block<3 || h->data_type<1 || h->data_type>7) return(1);
00347   if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
00348     tovax=1; else tovax=0;
00349 
00350   /* Copy short ints to buf */
00351   memcpy(buf+126, &h->data_type, 2);
00352   memcpy(buf+132, &h->dimension_1, 2); memcpy(buf+134, &h->dimension_2, 2);
00353   memcpy(buf+372, &h->norm_hour, 2);
00354   memcpy(buf+376, &h->norm_minute, 2);
00355   memcpy(buf+380, &h->norm_second, 2);
00356   memcpy(buf+384, &h->norm_day, 2);
00357   memcpy(buf+388, &h->norm_month, 2);
00358   memcpy(buf+392, &h->norm_year, 2);
00359   /* big to little endian if necessary */
00360   if(!little) swabip(buf, MatBLKSIZE);
00361 
00362   /* Copy floats to buf */
00363   ecat63wFloat(&h->scale_factor, buf+182, tovax, little);
00364   ecat63wFloat(&h->fov_source_width, buf+198, tovax, little);
00365 
00366   /* Write subheader */
00367   fseek(fp, (block-1)*MatBLKSIZE, SEEK_SET);
00368   if(ftell(fp)!=(block-1)*MatBLKSIZE) return(2);
00369   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
00370 
00371   return(0);
00372 }
00373 /*****************************************************************************/
00374 
00375 /*****************************************************************************/
00386 FILE *ecat63Create(const char *fname, ECAT63_mainheader *h) {
00387   FILE *fp;
00388   char tmp[FILENAME_MAX];
00389   int buf[MatBLKSIZE/4];
00390 
00391   if(ECAT63_TEST) printf("ecat63Create()\n");
00392   /* Check the arguments */
00393   if(fname==NULL || h==NULL) return(NULL);
00394   /* Check if file exists; backup, if necessary */
00395   if(access(fname, 0) != -1) {
00396     strcpy(tmp, fname); strcat(tmp, BACKUP_EXTENSION);
00397     if(access(tmp, 0) != -1) remove(tmp);
00398     if(ECAT63_TEST) printf("Renaming %s -> %s\n", fname, tmp);
00399     rename(fname, tmp);
00400   }
00401   /* Open file */
00402   fp=fopen(fname, "wb+"); if(fp==NULL) return(fp);
00403   /* Write main header */
00404   if(ecat63WriteMainheader(fp, h)) return(NULL);
00405   /* Construct an empty matrix list ; convert to little endian if necessary */
00406   memset(buf, 0, MatBLKSIZE);  
00407   buf[0]=31; buf[1]=2; if(!little_endian()) swawbip(buf, MatBLKSIZE);
00408   /* Write data buffer */
00409   fseek(fp, (MatFirstDirBlk-1)*MatBLKSIZE, SEEK_SET);
00410   if(ftell(fp)!=(MatFirstDirBlk-1)*MatBLKSIZE) return(NULL);
00411   if(fwrite(buf, 4, MatBLKSIZE/4, fp) != MatBLKSIZE/4) return(NULL);
00412   /* OK, then return file pointer */
00413   return(fp);
00414 }
00415 /*****************************************************************************/
00416 
00417 /*****************************************************************************/
00429 int ecat63WriteImage(FILE *fp, int matnum, ECAT63_imageheader *h, void *data) {
00430   int nxtblk, blkNr, data_size, pxlNr, pxlSize, ret;
00431 
00432   if(ECAT63_TEST) printf("ecat63WriteImage(fp, %d, ih, data)\n", matnum);
00433   if(fp==NULL || matnum<1 || h==NULL || data==NULL) return(1);
00434   /* nr of pixels */
00435   pxlNr=h->dimension_1*h->dimension_2; if(pxlNr<1) return(2);
00436   /* mem taken by one pixel */
00437   switch(h->data_type) {
00438     case BYTE_TYPE: pxlSize=1;
00439       break;
00440     case VAX_I2:
00441     case SUN_I2: pxlSize=2;
00442       break;
00443     case VAX_I4: return(3);
00444     case VAX_R4: return(3);
00445     case IEEE_R4:
00446     case SUN_I4: pxlSize=4;
00447       break;
00448     default: return(2);
00449   }
00450   /* mem taken by all pixels */
00451   data_size=pxlNr*pxlSize;
00452   /* block nr taken by all pixels */
00453   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(3);
00454   /* Get block number for matrix header and data */
00455   nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) return(4);
00456   if(ECAT63_TEST) printf("  block=%d\n", nxtblk);
00457   /* Write header */
00458   ret=ecat63WriteImageheader(fp, nxtblk, h); if(ret) return(40+ret);
00459   /* Write matrix data */
00460   ret=ecat63WriteMatdata(fp, nxtblk+1, data, pxlNr, pxlSize);
00461   if(ret) return(50+ret);
00462   return 0;
00463 }
00464 /*****************************************************************************/
00465 
00466 /*****************************************************************************/
00478 int ecat63WriteScan(FILE *fp, int matnum, ECAT63_scanheader *h, void *data) {
00479   int nxtblk, blkNr, data_size, pxlNr, pxlSize, ret;
00480 
00481   if(ECAT63_TEST) printf("ecat63WriteScan(fp, %d, sh, data)\n", matnum);
00482   if(fp==NULL || matnum<1 || h==NULL || data==NULL) return(1);
00483   /* nr of pixels */
00484   pxlNr=h->dimension_1*h->dimension_2; if(pxlNr<1) return(1);
00485   /* mem taken by one pixel */
00486   switch(h->data_type) {
00487     case BYTE_TYPE: pxlSize=1;
00488       break;
00489     case VAX_I2:
00490     case SUN_I2: pxlSize=2;
00491       break;
00492     case VAX_I4: return(3);
00493     case VAX_R4: return(3);
00494     case IEEE_R4:
00495     case SUN_I4: pxlSize=4;
00496       break;
00497     default: return(2);
00498   }
00499   /* mem taken by all pixels */
00500   data_size=pxlNr*pxlSize;
00501   /* block nr taken by all pixels */
00502   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(3);
00503   /* Get block number for matrix header and data */
00504   nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) return(4);
00505   if(ECAT63_TEST) printf("  block=%d\n", nxtblk);
00506   /* Write header */
00507   ret=ecat63WriteScanheader(fp, nxtblk, h); if(ret) return(40+ret);
00508   /* Write matrix data */
00509   ret=ecat63WriteMatdata(fp, nxtblk+1, data, pxlNr, pxlSize);
00510   if(ret) return(50+ret);
00511   return 0;
00512 }
00513 /*****************************************************************************/
00514 
00515 /*****************************************************************************/
00527 int ecat63WriteNorm(FILE *fp, int matnum, ECAT63_normheader *h, void *data) {
00528   int nxtblk, blkNr, data_size, pxlNr, pxlSize, ret;
00529 
00530   if(ECAT63_TEST) printf("ecat63WriteNorm(fp, %d, nh, data)\n", matnum);
00531   if(fp==NULL || matnum<1 || h==NULL || data==NULL) return(1);
00532   /* nr of pixels */
00533   pxlNr=h->dimension_1*h->dimension_2; if(pxlNr<1) return(1);
00534   /* mem taken by one pixel */
00535   switch(h->data_type) {
00536     case BYTE_TYPE: pxlSize=1;
00537       break;
00538     case VAX_I2:
00539     case SUN_I2: pxlSize=2;
00540       break;
00541     case VAX_I4: return(3);
00542     case VAX_R4: return(3);
00543     case IEEE_R4:
00544     case SUN_I4: pxlSize=4;
00545       break;
00546     default: return(2);
00547   }
00548   /* mem taken by all pixels */
00549   data_size=pxlNr*pxlSize;
00550   /* block nr taken by all pixels */
00551   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(3);
00552   /* Get block number for matrix header and data */
00553   nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) return(4);
00554   if(ECAT63_TEST) printf("  block=%d\n", nxtblk);
00555   /* Write header */
00556   ret=ecat63WriteNormheader(fp, nxtblk, h); if(ret) return(40+ret);
00557   /* Write matrix data */
00558   ret=ecat63WriteMatdata(fp, nxtblk+1, data, pxlNr, pxlSize);
00559   if(ret) return(50+ret);
00560   return 0;
00561 }
00562 /*****************************************************************************/
00563 
00564 /*****************************************************************************/
00576 int ecat63WriteAttn(FILE *fp, int matnum, ECAT63_attnheader *h, void *data) {
00577   int nxtblk, blkNr, data_size, pxlNr, pxlSize, ret;
00578 
00579   if(ECAT63_TEST) printf("ecat63WriteAttn(fp, %d, ah, data)\n", matnum);
00580   if(fp==NULL || matnum<1 || h==NULL || data==NULL) return(1);
00581   /* nr of pixels */
00582   pxlNr=h->dimension_1*h->dimension_2; if(pxlNr<1) return(1);
00583   /* mem taken by one pixel */
00584   switch(h->data_type) {
00585     case BYTE_TYPE: pxlSize=1;
00586       break;
00587     case VAX_I2:
00588     case SUN_I2: pxlSize=2;
00589       break;
00590     case VAX_I4: return(3);
00591     case VAX_R4: return(3);
00592     case IEEE_R4:
00593     case SUN_I4: pxlSize=4;
00594       break;
00595     default: return(2);
00596   }
00597   /* mem taken by all pixels */
00598   data_size=pxlNr*pxlSize;
00599   /* block nr taken by all pixels */
00600   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(3);
00601   /* Get block number for matrix header and data */
00602   nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) return(4);
00603   if(ECAT63_TEST) printf("  block=%d\n", nxtblk);
00604   /* Write header */
00605   ret=ecat63WriteAttnheader(fp, nxtblk, h); if(ret) return(40+ret);
00606   /* Write matrix data */
00607   ret=ecat63WriteMatdata(fp, nxtblk+1, data, pxlNr, pxlSize);
00608   if(ret) return(50+ret);
00609   return 0;
00610 }
00611 /*****************************************************************************/
00612 
00613 /*****************************************************************************/
00629 int ecat63WriteMatdata(FILE *fp, int strtblk, char *data, int pxlNr, int pxlSize) {
00630   unsigned char buf[MatBLKSIZE];
00631   char *dptr;
00632   int i, blkNr, dataSize, byteNr, little;
00633 
00634   if(ECAT63_TEST) printf("ecat63WriteMatdata(fp, %d, data, %d, %d)\n", strtblk, pxlNr, pxlSize);
00635   if(fp==NULL || strtblk<1 || data==NULL || pxlNr<1 || pxlSize<1) return(1);
00636   little=little_endian(); memset(buf, 0, MatBLKSIZE);
00637   dataSize=pxlNr*pxlSize; if(dataSize<1) return(1);
00638   /* block nr taken by all pixels */
00639   blkNr=(dataSize+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(1);
00640   if(ECAT63_TEST>1) printf("    blkNr=%d\n", blkNr);
00641   /* Search the place for writing */
00642   fseek(fp, (strtblk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(strtblk-1)*MatBLKSIZE) return(2);
00643   /* Save blocks one at a time */
00644   for(i=0, dptr=data; i<blkNr && dataSize>0; i++) {
00645     byteNr=(dataSize<MatBLKSIZE)?dataSize:MatBLKSIZE;
00646     memcpy(buf, dptr, byteNr);
00647     /* Change matrix byte order in big endian platforms */
00648     if(!little_endian()) {
00649       if(pxlSize==2) swabip(buf, byteNr);
00650       else if(pxlSize==4) swawbip(buf, byteNr);
00651     }
00652     /* Write block */
00653     if(fwrite(buf, 1, MatBLKSIZE, fp)!=MatBLKSIZE) return(3);
00654     /* Prepare for the next block */
00655     dptr+=byteNr;
00656     dataSize-=byteNr;
00657   } /* next block */
00658   return(0);
00659 }
00660 /*****************************************************************************/
00661 
00662 /*****************************************************************************/
00672 int ecat63_is_scaling_needed(float amax, float *data, int nr) {
00673   int i;
00674   double d;
00675 
00676   if(nr<1 || data==NULL) return(0);
00677   /* scaling is necessary if all values are between -1 - 1 */
00678   if(amax<0.9999) return(1);
00679   /* Lets check first if at least the max value is close to integers or not */
00680   if(modf(amax, &d)>0.0001) return(1);
00681   /* if it is, then check all pixels */
00682   for(i=0; i<nr; i++) if(modf(*data++, &d)>0.0001) return(1);
00683   return(0);
00684 }
00685 /*****************************************************************************/
00686 
00687 /*****************************************************************************/
00700 int ecat63WriteImageMatrix(FILE *fp, int matnum, ECAT63_imageheader *h, float *fdata) {
00701   int i, nxtblk, blkNr, data_size, pxlNr, ret;
00702   float *fptr, fmin, fmax, g, f;
00703   char *mdata, *mptr;
00704   short int *sptr;
00705 
00706 
00707 
00708   if(ECAT63_TEST) printf("ecat63WriteImageMatrix(fp, %d, h, data)\n", matnum);
00709   if(fp==NULL || matnum<1 || h==NULL || fdata==NULL) {
00710     sprintf(ecat63errmsg, "invalid function parameter.\n");
00711     return(1);
00712   }
00713   /* nr of pixels */
00714   pxlNr=h->dimension_1*h->dimension_2;
00715   if(pxlNr<1) {
00716     sprintf(ecat63errmsg, "invalid matrix dimension.\n");
00717     return(3);
00718   }
00719   /* How much memory is needed for ALL pixels */
00720   data_size=pxlNr*ecat63pxlbytes(h->data_type);
00721   /* block nr taken by all pixels */
00722   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) {
00723     sprintf(ecat63errmsg, "invalid block number.\n");
00724     return(4);
00725   }
00726   /* Allocate memory for matrix data */
00727   mdata=(char*)calloc(blkNr, MatBLKSIZE); if(mdata==NULL) {
00728     sprintf(ecat63errmsg, "out of memory.\n");
00729     return(5);
00730   }
00731   /* Search for min and max for calculation of scale factor */
00732   fptr=fdata; fmin=fmax=*fptr;
00733   for(i=0; i<pxlNr; i++, fptr++) {
00734     if(*fptr>fmax) fmax=*fptr; else if(*fptr<fmin) fmin=*fptr;
00735   }
00736   if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
00737   if(g>0) f=32766./g; else f=1.0;
00738   /* Check if pixels values can be left as such with scale_factor = 1 */
00739   fptr=fdata;
00740   if(f>=1.0 && ecat63_is_scaling_needed(g, fptr, pxlNr)==0) f=1.0;
00741   /* Scale matrix data to shorts */
00742   h->quant_scale=1.0/f;
00743   sptr=(short int*)mdata; fptr=fdata;
00744   for(i=0; i<pxlNr; i++, sptr++, fptr++) *sptr=(short int)temp_roundf(f*(*fptr));
00745   /* Set header short min & max */
00746   h->image_min=(short int)temp_roundf(f*fmin);
00747   h->image_max=(short int)temp_roundf(f*fmax);
00748   /* Get block number for matrix header and data */
00749   nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) {
00750     sprintf(ecat63errmsg, "cannot determine matrix block (%d).\n", -nxtblk);
00751     free(mdata); return(8);
00752   }
00753   if(ECAT63_TEST>2) printf("  block=%d fmin=%g fmax=%g\n", nxtblk, fmin, fmax);
00754   /* Write header */
00755   ret=ecat63WriteImageheader(fp, nxtblk, h); if(ret) {
00756     sprintf(ecat63errmsg, "cannot write subheader (%d).\n", ret);
00757     free(mdata); return(10);
00758   }
00759   /* Write matrix data */
00760   mptr=mdata;
00761   ret=ecat63WriteMatdata(fp, nxtblk+1, mptr, pxlNr, ecat63pxlbytes(h->data_type));
00762   free(mdata);
00763   if(ret) {
00764     sprintf(ecat63errmsg, "cannot write matrix data (%d).\n", ret);
00765     return(13);
00766   }
00767   return(0);
00768 }
00769 /*****************************************************************************/
00770 
00771 /*****************************************************************************/
00784 int ecat63WriteScanMatrix(FILE *fp, int matnum, ECAT63_scanheader *h, float *fdata) {
00785   int i, nxtblk, blkNr, data_size, pxlNr, ret;
00786   float *fptr, fmin, fmax, g, f;
00787   char *mdata, *mptr;
00788   short int *sptr;
00789 
00790 
00791   if(ECAT63_TEST) printf("ecat63WriteScanMatrix(fp, %d, h, data)\n", matnum);
00792   if(fp==NULL || matnum<1 || h==NULL || fdata==NULL) {
00793     sprintf(ecat63errmsg, "invalid function parameter.\n");
00794     return(1);
00795   }
00796   /* nr of pixels */
00797   pxlNr=h->dimension_1*h->dimension_2;
00798   if(pxlNr<1) {
00799     sprintf(ecat63errmsg, "invalid matrix dimension.\n");
00800     return(3);
00801   }
00802   /* How much memory is needed for ALL pixels */
00803   data_size=pxlNr*ecat63pxlbytes(h->data_type);
00804   /* block nr taken by all pixels */
00805   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) {
00806     sprintf(ecat63errmsg, "invalid block number.\n");
00807     return(4);
00808   }
00809   /* Allocate memory for matrix data */
00810   mdata=(char*)calloc(blkNr, MatBLKSIZE); if(mdata==NULL) {
00811     sprintf(ecat63errmsg, "out of memory.\n");
00812     return(5);
00813   }
00814   /* Search for min and max for calculation of scale factor */
00815   fptr=fdata; fmin=fmax=*fptr;
00816   for(i=0; i<pxlNr; i++, fptr++) {
00817     if(*fptr>fmax) fmax=*fptr; else if(*fptr<fmin) fmin=*fptr;
00818   }
00819   if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
00820   if(g>0) f=32766./g; else f=1.0;
00821   /* Check if pixels values can be left as such with scale_factor = 1 */
00822   fptr=fdata;
00823   if(f>=1.0 && ecat63_is_scaling_needed(g, fptr, pxlNr)==0) f=1.0;
00824   /* Scale matrix data to shorts */
00825   h->scale_factor=1.0/f;
00826   sptr=(short int*)mdata; fptr=fdata;
00827   for(i=0; i<pxlNr; i++, sptr++, fptr++) *sptr=(short int)temp_roundf(f*(*fptr));
00828   /* Set header short min & max */
00829   h->scan_min=(short int)temp_roundf(f*fmin);
00830   h->scan_max=(short int)temp_roundf(f*fmax);
00831   /* Get block number for matrix header and data */
00832   nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) {
00833     sprintf(ecat63errmsg, "cannot determine matrix block (%d).\n", -nxtblk);
00834     free(mdata); return(8);
00835   }
00836   if(ECAT63_TEST>2) printf("  block=%d fmin=%g fmax=%g\n", nxtblk, fmin, fmax);
00837   /* Write header */
00838   ret=ecat63WriteScanheader(fp, nxtblk, h); if(ret) {
00839     sprintf(ecat63errmsg, "cannot write subheader (%d).\n", ret);
00840     free(mdata); return(10);
00841   }
00842   /* Write matrix data */
00843   mptr=mdata;
00844   ret=ecat63WriteMatdata(fp, nxtblk+1, mptr, pxlNr, ecat63pxlbytes(h->data_type));
00845   free(mdata);
00846   if(ret) {
00847     sprintf(ecat63errmsg, "cannot write matrix data (%d).\n", ret);
00848     return(13);
00849   }
00850   return(0);
00851 }
00852 /*****************************************************************************/
00853 
00854 /*****************************************************************************/
00863 void ecat63wFloat(float *bufi, void *bufo, int tovax, int islittle) {
00864   unsigned int ul;
00865 
00866   memcpy(&ul, bufi, 4); if(ul==0) {memcpy(bufo, bufi, 4); return;}
00867   if(tovax) { /* If VAX format is needed */
00868     ul+=(2L<<23); /* increase exp by 2 */
00869     /* Swap words on i386 and bytes on SUN */
00870     if(islittle) swawip(&ul, 4); else swabip(&ul, 4);
00871   } else {
00872     if(!islittle) swawbip(&ul, 4); /* Switch words and bytes on SUN */
00873   }
00874   memcpy(bufo, &ul, 4);
00875 }
00885 void ecat63wInt(int *bufi, void *bufo, int tovax, int islittle) {
00886   int i;
00887 
00888   /* Swap both words and bytes on SUN */
00889   memcpy(&i, bufi, 4); if(!islittle) swawbip(&i, 4);
00890   memcpy(bufo, &i, 4);
00891 }
00892 /*****************************************************************************/
00893 
00894 /*****************************************************************************/