ecat7w.c
Go to the documentation of this file.
00001 /******************************************************************************
00002 
00003   Copyright (c) 2003-2008 Turku PET Centre
00004 
00005   Library file: ecat7w.c
00006   Description:  Functions for writing ECAT 7.x format.
00007 
00008   This program is free software; you can redistribute it and/or modify it under
00009   the terms of the GNU General Public License as published by the Free Software
00010   Foundation; either version 2 of the License, or (at your option) any later
00011   version.
00012 
00013   This program is distributed in the hope that it will be useful, but WITHOUT
00014   ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
00015   FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
00016 
00017   You should have received a copy of the GNU General Public License along with
00018   this program; if not, write to the Free Software Foundation, Inc., 59 Temple
00019   Place, Suite 330, Boston, MA 02111-1307 USA.
00020 
00021   Turku PET Centre hereby disclaims all copyright interest in the program.
00022   Juhani Knuuti
00023   Director, Professor
00024   Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
00025   
00026   Modification history:
00027   2003-07-24 Vesa Oikonen
00028     First created.
00029   2003-09-04 VO
00030     Added support for 3D sinograms, ecat7WriteScanMatrix().
00031     Removed functions ecat7wFloat() and ecat7wInt().
00032   2003-11-30 VO
00033     For now, calls temp_roundf() instead of roundf().
00034   2004-01-07 VO
00035     ecat7WriteImageMatrix(): corrected img min&max in header.
00036   2004-09-17 VO
00037     Changes in comment style.
00038     Removed code that was previously commented out.
00039   2004-12-28 VO
00040     Included function ecat7_is_scaling_needed().
00041     This function is applied to determine whether scal factor can be set to
00042     one in case that all pixel values are close to integers and small enough.
00043   2007-01-27 VO
00044     Unsigned char pointer was corrected to signed in ecat7WriteMatrixdata().
00045   2007-03-27 VO
00046     Added ecat7WritePolarmapMatrix().
00047   2007-09-02 VO
00048     Backup file extension changed from % to .bak.
00049 
00050 ******************************************************************************/
00051 #include <stdio.h>
00052 #include <stdlib.h>
00053 #include <math.h>
00054 #include <string.h>
00055 #include <unistd.h>
00056 #include <time.h>
00057 /*****************************************************************************/
00058 #include "swap.h"
00059 #include "petc99.h"
00060 #include "include/ecat7.h"
00061 /*****************************************************************************/
00062 
00063 /*****************************************************************************/
00073 int ecat7WriteMainheader(FILE *fp, ECAT7_mainheader *h) {
00074   unsigned char buf[MatBLKSIZE];
00075   int little;
00076 
00077   if(ECAT7_TEST) printf("ecat7WriteMainheader()\n");
00078   /* Check arguments */
00079   if(fp==NULL || h==NULL) return(1);
00080   little=little_endian();
00081   /* Clear buf */
00082   memset(buf, 0, MatBLKSIZE);
00083 
00084   /* Copy header contents into buffer and change byte order if necessary */
00085   memcpy(buf+0, &h->magic_number, 14);
00086   memcpy(buf+14, &h->original_file_name, 32);
00087   memcpy(buf+46, &h->sw_version, 2); if(little) swabip(buf+46, 2);
00088   memcpy(buf+48, &h->system_type, 2); if(little) swabip(buf+48, 2);
00089   memcpy(buf+50, &h->file_type, 2); if(little) swabip(buf+50, 2);
00090   memcpy(buf+52, &h->serial_number, 10);
00091   memcpy(buf+62, &h->scan_start_time, 4); if(little) swawbip(buf+62, 4);
00092   memcpy(buf+66, &h->isotope_name, 8);
00093   memcpy(buf+74, &h->isotope_halflife, 4); if(little) swawbip(buf+74, 4);
00094   memcpy(buf+78, &h->radiopharmaceutical, 32);
00095   memcpy(buf+110, &h->gantry_tilt, 4); if(little) swawbip(buf+110, 4);
00096   memcpy(buf+114, &h->gantry_rotation, 4); if(little) swawbip(buf+114, 4);
00097   memcpy(buf+118, &h->bed_elevation, 4); if(little) swawbip(buf+118, 4);
00098   memcpy(buf+122, &h->intrinsic_tilt, 4); if(little) swawbip(buf+122, 4);
00099   memcpy(buf+126, &h->wobble_speed, 2); if(little) swabip(buf+126, 2);
00100   memcpy(buf+128, &h->transm_source_type, 2); if(little) swabip(buf+128, 2);
00101   memcpy(buf+130, &h->distance_scanned, 4); if(little) swawbip(buf+130, 4);
00102   memcpy(buf+134, &h->transaxial_fov, 4); if(little) swawbip(buf+134, 4);
00103   memcpy(buf+138, &h->angular_compression, 2); if(little) swabip(buf+138, 2);
00104   memcpy(buf+140, &h->coin_samp_mode, 2); if(little) swabip(buf+140, 2);
00105   memcpy(buf+142, &h->axial_samp_mode, 2); if(little) swabip(buf+142, 2);
00106   memcpy(buf+144, &h->ecat_calibration_factor, 4); if(little) swawbip(buf+144, 4);
00107   memcpy(buf+148, &h->calibration_units, 2); if(little) swabip(buf+148, 2);
00108   memcpy(buf+150, &h->calibration_units_label, 2); if(little) swabip(buf+150, 2);
00109   memcpy(buf+152, &h->compression_code, 2); if(little) swabip(buf+152, 2);
00110   memcpy(buf+154, &h->study_type, 12);
00111   memcpy(buf+166, &h->patient_id, 16);
00112   memcpy(buf+182, &h->patient_name, 32);
00113   memcpy(buf+214, &h->patient_sex, 1);
00114   memcpy(buf+215, &h->patient_dexterity, 1);
00115   memcpy(buf+216, &h->patient_age, 4); if(little) swawbip(buf+216, 4);
00116   memcpy(buf+220, &h->patient_height, 4); if(little) swawbip(buf+220, 4);
00117   memcpy(buf+224, &h->patient_weight, 4); if(little) swawbip(buf+224, 4);
00118   memcpy(buf+228, &h->patient_birth_date, 4); if(little) swawbip(buf+228, 4);
00119   memcpy(buf+232, &h->physician_name, 32);
00120   memcpy(buf+264, &h->operator_name, 32);
00121   memcpy(buf+296, &h->study_description, 32);
00122   memcpy(buf+328, &h->acquisition_type, 2); if(little) swabip(buf+328, 2);
00123   memcpy(buf+330, &h->patient_orientation, 2); if(little) swabip(buf+330, 2);
00124   memcpy(buf+332, &h->facility_name, 20);
00125   memcpy(buf+352, &h->num_planes, 2); if(little) swabip(buf+352, 2);
00126   memcpy(buf+354, &h->num_frames, 2); if(little) swabip(buf+354, 2);
00127   memcpy(buf+356, &h->num_gates, 2); if(little) swabip(buf+356, 2);
00128   memcpy(buf+358, &h->num_bed_pos, 2); if(little) swabip(buf+358, 2);
00129   memcpy(buf+360, &h->init_bed_position, 4); if(little) swawbip(buf+360, 4);
00130   memcpy(buf+364, h->bed_position, 15*4); if(little) swawbip(buf+364, 15*4);
00131   memcpy(buf+424, &h->plane_separation, 4); if(little) swawbip(buf+424, 4);
00132   memcpy(buf+428, &h->lwr_sctr_thres, 2); if(little) swabip(buf+428, 2);
00133   memcpy(buf+430, &h->lwr_true_thres, 2); if(little) swabip(buf+430, 2);
00134   memcpy(buf+432, &h->upr_true_thres, 2); if(little) swabip(buf+432, 2);
00135   memcpy(buf+434, &h->user_process_code, 10);
00136   memcpy(buf+444, &h->acquisition_mode, 2); if(little) swabip(buf+444, 2);
00137   memcpy(buf+446, &h->bin_size, 4); if(little) swawbip(buf+446, 4);
00138   memcpy(buf+450, &h->branching_fraction, 4); if(little) swawbip(buf+450, 4);
00139   memcpy(buf+454, &h->dose_start_time, 4); if(little) swawbip(buf+454, 4);
00140   memcpy(buf+458, &h->dosage, 4); if(little) swawbip(buf+458, 4);
00141   memcpy(buf+462, &h->well_counter_corr_factor, 4); if(little) swawbip(buf+462, 4);
00142   memcpy(buf+466, &h->data_units, 32);
00143   memcpy(buf+498, &h->septa_state, 2); if(little) swabip(buf+498, 2);
00144   memcpy(buf+500, &h->fill_cti, 12);
00145 
00146   /* Write main header */
00147   fseek(fp, 0*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=0*MatBLKSIZE) return(4);
00148   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(5);
00149 
00150   return(0);
00151 }
00152 /*****************************************************************************/
00153 
00154 /*****************************************************************************/
00164 int ecat7WriteImageheader(FILE *fp, int blk, ECAT7_imageheader *h) {
00165   unsigned char buf[MatBLKSIZE];
00166   int little; /* 1 if current platform is little endian (i386), else 0 */
00167 
00168   if(ECAT7_TEST) printf("ecat7WriteImageheader()\n");
00169   if(fp==NULL || blk<2 || h==NULL) return(1);
00170   little=little_endian(); if(ECAT7_TEST) printf("little=%d\n", little);
00171   /* Clear buf */
00172   memset(buf, 0, MatBLKSIZE);
00173   if(h->data_type==ECAT7_VAXI2) h->data_type=ECAT7_SUNI2;
00174   else if(h->data_type==ECAT7_VAXI4) h->data_type=ECAT7_SUNI4;
00175   else if(h->data_type==ECAT7_VAXR4) h->data_type=ECAT7_IEEER4;
00176   
00177   /* Copy the header fields and swap if necessary */
00178   memcpy(buf+0, &h->data_type, 2); if(little) swabip(buf+0, 2);
00179   memcpy(buf+2, &h->num_dimensions, 2); if(little) swabip(buf+2, 2);
00180   memcpy(buf+4, &h->x_dimension, 2); if(little) swabip(buf+4, 2);
00181   memcpy(buf+6, &h->y_dimension, 2); if(little) swabip(buf+6, 2);
00182   memcpy(buf+8, &h->z_dimension, 2); if(little) swabip(buf+8, 2);
00183   memcpy(buf+10, &h->x_offset, 4); if(little) swawbip(buf+10, 4);
00184   memcpy(buf+14, &h->y_offset, 4); if(little) swawbip(buf+14, 4);
00185   memcpy(buf+18, &h->z_offset, 4); if(little) swawbip(buf+18, 4);
00186   memcpy(buf+22, &h->recon_zoom, 4); if(little) swawbip(buf+22, 4);
00187   memcpy(buf+26, &h->scale_factor, 4); if(little) swawbip(buf+26, 4);
00188   memcpy(buf+30, &h->image_min, 2); if(little) swabip(buf+30, 2);
00189   memcpy(buf+32, &h->image_max, 2); if(little) swabip(buf+32, 2);
00190   memcpy(buf+34, &h->x_pixel_size, 4); if(little) swawbip(buf+34, 4);
00191   memcpy(buf+38, &h->y_pixel_size, 4); if(little) swawbip(buf+38, 4);
00192   memcpy(buf+42, &h->z_pixel_size, 4); if(little) swawbip(buf+42, 4);
00193   memcpy(buf+46, &h->frame_duration, 4); if(little) swawbip(buf+46, 4);
00194   memcpy(buf+50, &h->frame_start_time, 4); if(little) swawbip(buf+50, 4);
00195   memcpy(buf+54, &h->filter_code, 2); if(little) swabip(buf+54, 2);
00196   memcpy(buf+56, &h->x_resolution, 4); if(little) swawbip(buf+56, 4);
00197   memcpy(buf+60, &h->y_resolution, 4); if(little) swawbip(buf+60, 4);
00198   memcpy(buf+64, &h->z_resolution, 4); if(little) swawbip(buf+64, 4);
00199   memcpy(buf+68, &h->num_r_elements, 4); if(little) swawbip(buf+68, 4);
00200   memcpy(buf+72, &h->num_angles, 4); if(little) swawbip(buf+72, 4);
00201   memcpy(buf+76, &h->z_rotation_angle, 4); if(little) swawbip(buf+76, 4);
00202   memcpy(buf+80, &h->decay_corr_fctr, 4); if(little) swawbip(buf+80, 4);
00203   memcpy(buf+84, &h->processing_code, 4); if(little) swawbip(buf+84, 4);
00204   memcpy(buf+88, &h->gate_duration, 4); if(little) swawbip(buf+88, 4);
00205   memcpy(buf+92, &h->r_wave_offset, 4); if(little) swawbip(buf+92, 4);
00206   memcpy(buf+96, &h->num_accepted_beats, 4); if(little) swawbip(buf+96, 4);
00207   memcpy(buf+100, &h->filter_cutoff_frequency, 4); if(little) swawbip(buf+100, 4);
00208   memcpy(buf+104, &h->filter_resolution, 4); if(little) swawbip(buf+104, 4);
00209   memcpy(buf+108, &h->filter_ramp_slope, 4); if(little) swawbip(buf+108, 4);
00210   memcpy(buf+112, &h->filter_order, 2); if(little) swabip(buf+112, 2);
00211   memcpy(buf+114, &h->filter_scatter_fraction, 4); if(little) swawbip(buf+114, 4);
00212   memcpy(buf+118, &h->filter_scatter_slope, 4); if(little) swawbip(buf+118, 4);
00213   memcpy(buf+122, &h->annotation, 40);
00214   memcpy(buf+162, &h->mt_1_1, 4); if(little) swawbip(buf+162, 4);
00215   memcpy(buf+166, &h->mt_1_2, 4); if(little) swawbip(buf+166, 4);
00216   memcpy(buf+170, &h->mt_1_3, 4); if(little) swawbip(buf+170, 4);
00217   memcpy(buf+174, &h->mt_2_1, 4); if(little) swawbip(buf+174, 4);
00218   memcpy(buf+178, &h->mt_2_2, 4); if(little) swawbip(buf+178, 4);
00219   memcpy(buf+182, &h->mt_2_3, 4); if(little) swawbip(buf+182, 4);
00220   memcpy(buf+186, &h->mt_3_1, 4); if(little) swawbip(buf+186, 4);
00221   memcpy(buf+190, &h->mt_3_2, 4); if(little) swawbip(buf+190, 4);
00222   memcpy(buf+194, &h->mt_3_3, 4); if(little) swawbip(buf+194, 4);
00223   memcpy(buf+198, &h->rfilter_cutoff, 4); if(little) swawbip(buf+198, 4);
00224   memcpy(buf+202, &h->rfilter_resolution, 4); if(little) swawbip(buf+202, 4);
00225   memcpy(buf+206, &h->rfilter_code, 2); if(little) swabip(buf+206, 2);
00226   memcpy(buf+208, &h->rfilter_order, 2); if(little) swabip(buf+208, 2);
00227   memcpy(buf+210, &h->zfilter_cutoff, 4); if(little) swawbip(buf+210, 4);
00228   memcpy(buf+214, &h->zfilter_resolution, 4); if(little) swawbip(buf+214, 4);
00229   memcpy(buf+218, &h->zfilter_code, 2); if(little) swabip(buf+218, 2);
00230   memcpy(buf+220, &h->zfilter_order, 2); if(little) swabip(buf+220, 2);
00231   memcpy(buf+222, &h->mt_1_4, 4); if(little) swawbip(buf+222, 4);
00232   memcpy(buf+226, &h->mt_2_4, 4); if(little) swawbip(buf+226, 4);
00233   memcpy(buf+230, &h->mt_3_4, 4); if(little) swawbip(buf+230, 4);
00234   memcpy(buf+234, &h->scatter_type, 2); if(little) swabip(buf+234, 2);
00235   memcpy(buf+236, &h->recon_type, 2); if(little) swabip(buf+236, 2);
00236   memcpy(buf+238, &h->recon_views, 2); if(little) swabip(buf+238, 2);
00237   memcpy(buf+240, &h->fill_cti, 87);
00238   memcpy(buf+414, &h->fill_user, 48);
00239 
00240   /* Write header */
00241   fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(4);
00242   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(5);
00243 
00244   return(0);
00245 }
00246 /*****************************************************************************/
00247 
00248 /*****************************************************************************/
00258 int ecat7WriteAttenheader(FILE *fp, int blk, ECAT7_attenheader *h) {
00259   unsigned char buf[MatBLKSIZE];
00260   int little; /* 1 if current platform is little endian (i386), else 0 */
00261 
00262   if(ECAT7_TEST) printf("ecat7WriteAttenheader()\n");
00263   if(fp==NULL || blk<2 || h==NULL) return(1);
00264   little=little_endian(); if(ECAT7_TEST) printf("little=%d\n", little);
00265   /* Clear buf */
00266   memset(buf, 0, MatBLKSIZE);
00267   if(h->data_type==ECAT7_VAXI2) h->data_type=ECAT7_SUNI2;
00268   else if(h->data_type==ECAT7_VAXI4) h->data_type=ECAT7_SUNI4;
00269   else if(h->data_type==ECAT7_VAXR4) h->data_type=ECAT7_IEEER4;
00270   
00271   /* Copy the header fields and swap if necessary */
00272   memcpy(buf+0, &h->data_type, 2); if(little) swabip(buf+0, 2);
00273   memcpy(buf+2, &h->num_dimensions, 2); if(little) swabip(buf+2, 2);
00274   memcpy(buf+4, &h->attenuation_type, 2); if(little) swabip(buf+4, 2);
00275   memcpy(buf+6, &h->num_r_elements, 2); if(little) swabip(buf+6, 2);
00276   memcpy(buf+8, &h->num_angles, 2); if(little) swabip(buf+8, 2);
00277   memcpy(buf+10, &h->num_z_elements, 2); if(little) swabip(buf+10, 2);
00278   memcpy(buf+12, &h->ring_difference, 2); if(little) swabip(buf+12, 2);
00279   memcpy(buf+14, &h->x_resolution, 4); if(little) swawbip(buf+14, 4);
00280   memcpy(buf+18, &h->y_resolution, 4); if(little) swawbip(buf+18, 4);
00281   memcpy(buf+22, &h->z_resolution, 4); if(little) swawbip(buf+22, 4);
00282   memcpy(buf+26, &h->w_resolution, 4); if(little) swawbip(buf+26, 4);
00283   memcpy(buf+30, &h->scale_factor, 4); if(little) swawbip(buf+30, 4);
00284   memcpy(buf+34, &h->x_offset, 4); if(little) swawbip(buf+34, 4);
00285   memcpy(buf+38, &h->y_offset, 4); if(little) swawbip(buf+38, 4);
00286   memcpy(buf+42, &h->x_radius, 4); if(little) swawbip(buf+42, 4);
00287   memcpy(buf+46, &h->y_radius, 4); if(little) swawbip(buf+46, 4);
00288   memcpy(buf+50, &h->tilt_angle, 4); if(little) swawbip(buf+50, 4);
00289   memcpy(buf+54, &h->attenuation_coeff, 4); if(little) swawbip(buf+54, 4);
00290   memcpy(buf+58, &h->attenuation_min, 4); if(little) swawbip(buf+58, 4);
00291   memcpy(buf+62, &h->attenuation_max, 4); if(little) swawbip(buf+62, 4);
00292   memcpy(buf+66, &h->skull_thickness, 4); if(little) swawbip(buf+66, 4);
00293   memcpy(buf+70, &h->num_additional_atten_coeff, 2); if(little) swabip(buf+70, 2);
00294   memcpy(buf+72, h->additional_atten_coeff, 8*4); if(little) swawbip(buf+72, 8*4);
00295   memcpy(buf+104, &h->edge_finding_threshold, 4); if(little) swawbip(buf+104, 4);
00296   memcpy(buf+108, &h->storage_order, 2); if(little) swabip(buf+108, 2);
00297   memcpy(buf+110, &h->span, 2); if(little) swabip(buf+110, 2);
00298   memcpy(buf+112, h->z_elements, 64*2); if(little) swabip(buf+112, 64*2);
00299   memcpy(buf+240, h->fill_cti, 86*2); if(little) swabip(buf+240, 86*2);
00300   memcpy(buf+412, h->fill_user, 50*2); if(little) swabip(buf+412, 50*2);
00301 
00302   /* Write header */
00303   fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET);
00304   if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(4);
00305   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(5);
00306 
00307   return(0);
00308 }
00309 /*****************************************************************************/
00310 
00311 /*****************************************************************************/
00321 int ecat7WritePolmapheader(FILE *fp, int blk, ECAT7_polmapheader *h) {
00322   unsigned char buf[MatBLKSIZE];
00323   int little; /* 1 if current platform is little endian (i386), else 0 */
00324 
00325   if(ECAT7_TEST) printf("ecat7WritePolmapheader()\n");
00326   if(fp==NULL || blk<2 || h==NULL) return(1);
00327   little=little_endian(); if(ECAT7_TEST) printf("little=%d\n", little);
00328   /* Clear buf */
00329   memset(buf, 0, MatBLKSIZE);
00330   if(h->data_type==ECAT7_VAXI2) h->data_type=ECAT7_SUNI2;
00331   else if(h->data_type==ECAT7_VAXI4) h->data_type=ECAT7_SUNI4;
00332   else if(h->data_type==ECAT7_VAXR4) h->data_type=ECAT7_IEEER4;
00333   
00334   /* Copy the header fields and swap if necessary */
00335   memcpy(buf+0, &h->data_type, 2); if(little) swabip(buf+0, 2);
00336   memcpy(buf+2, &h->polar_map_type, 2); if(little) swabip(buf+2, 2);
00337   memcpy(buf+4, &h->num_rings, 2); if(little) swabip(buf+4, 2);
00338   memcpy(buf+6, h->sectors_per_ring, 32*2); if(little) swabip(buf+6, 32*2);
00339   memcpy(buf+70, h->ring_position, 32*4); if(little) swawbip(buf+70, 32*4);
00340   memcpy(buf+198, h->ring_angle, 32*2); if(little) swabip(buf+198, 32*2);
00341   memcpy(buf+262, &h->start_angle, 2); if(little) swabip(buf+262, 2);
00342   memcpy(buf+264, h->long_axis_left, 3*2); if(little) swabip(buf+264, 3*2);
00343   memcpy(buf+270, h->long_axis_right, 3*2); if(little) swabip(buf+270, 3*2);
00344   memcpy(buf+276, &h->position_data, 2); if(little) swabip(buf+276, 2);
00345   memcpy(buf+278, &h->image_min, 2); if(little) swabip(buf+278, 2);
00346   memcpy(buf+280, &h->image_max, 2); if(little) swabip(buf+280, 2);
00347   memcpy(buf+282, &h->scale_factor, 4); if(little) swawbip(buf+282, 4);
00348   memcpy(buf+286, &h->pixel_size, 4); if(little) swawbip(buf+286, 4);
00349   memcpy(buf+290, &h->frame_duration, 4); if(little) swawbip(buf+290, 4);
00350   memcpy(buf+294, &h->frame_start_time, 4); if(little) swawbip(buf+294, 4);
00351   memcpy(buf+298, &h->processing_code, 2); if(little) swabip(buf+298, 2);
00352   memcpy(buf+300, &h->quant_units, 2); if(little) swabip(buf+300, 2);
00353   memcpy(buf+302, h->annotation, 40);
00354   memcpy(buf+342, &h->gate_duration, 4); if(little) swawbip(buf+342, 4);
00355   memcpy(buf+346, &h->r_wave_offset, 4); if(little) swawbip(buf+346, 4);
00356   memcpy(buf+350, &h->num_accepted_beats, 4); if(little) swawbip(buf+350, 4);
00357   memcpy(buf+354, h->polar_map_protocol, 20);
00358   memcpy(buf+374, h->database_name, 30);
00359   memcpy(buf+404, h->fill_cti, 27*2); if(little) swabip(buf+404, 27*2);
00360 
00361   /* Write header */
00362   fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET);
00363   if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(4);
00364   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(5);
00365 
00366   return(0);
00367 }
00368 /*****************************************************************************/
00369 
00370 /*****************************************************************************/
00380 int ecat7WriteNormheader(FILE *fp, int blk, ECAT7_normheader *h) {
00381   unsigned char buf[MatBLKSIZE];
00382   int little; /* 1 if current platform is little endian (i386), else 0 */
00383 
00384   if(ECAT7_TEST) printf("ecat7WriteNormheader()\n");
00385   if(fp==NULL || blk<2 || h==NULL) return(1);
00386   little=little_endian(); if(ECAT7_TEST) printf("little=%d\n", little);
00387   /* Clear buf */
00388   memset(buf, 0, MatBLKSIZE);
00389   if(h->data_type==ECAT7_VAXI2) h->data_type=ECAT7_SUNI2;
00390   else if(h->data_type==ECAT7_VAXI4) h->data_type=ECAT7_SUNI4;
00391   else if(h->data_type==ECAT7_VAXR4) h->data_type=ECAT7_IEEER4;
00392   
00393   /* Copy the header fields and swap if necessary */
00394   memcpy(buf+0, &h->data_type, 2); if(little) swabip(buf+0, 2);
00395   memcpy(buf+2, &h->num_r_elements, 2); if(little) swabip(buf+2, 2);
00396   memcpy(buf+4, &h->num_transaxial_crystals, 2); if(little) swabip(buf+4, 2);
00397   memcpy(buf+6, &h->num_crystal_rings, 2); if(little) swabip(buf+6, 2);
00398   memcpy(buf+8, &h->crystals_per_ring, 2); if(little) swabip(buf+8, 2);
00399   memcpy(buf+10, &h->num_geo_corr_planes, 2); if(little) swabip(buf+10, 2);
00400   memcpy(buf+12, &h->uld, 2); if(little) swabip(buf+12, 2);
00401   memcpy(buf+14, &h->lld, 2); if(little) swabip(buf+14, 2);
00402   memcpy(buf+16, &h->scatter_energy, 2); if(little) swabip(buf+16, 2);
00403   memcpy(buf+18, &h->norm_quality_factor, 4); if(little) swawbip(buf+18, 4);
00404   memcpy(buf+22, &h->norm_quality_factor_code, 2); if(little) swabip(buf+22, 2);
00405   memcpy(buf+24, h->ring_dtcor1, 32*4); if(little) swawbip(buf+24, 32*4);
00406   memcpy(buf+152, h->ring_dtcor2, 32*4); if(little) swawbip(buf+152, 32*4);
00407   memcpy(buf+280, h->crystal_dtcor, 8*4); if(little) swawbip(buf+280, 8*4);
00408   memcpy(buf+312, &h->span, 2); if(little) swabip(buf+312, 2);
00409   memcpy(buf+314, &h->max_ring_diff, 2); if(little) swabip(buf+314, 2);
00410   memcpy(buf+316, h->fill_cti, 48*2); if(little) swabip(buf+316, 48*2);
00411   memcpy(buf+412, h->fill_user, 50*2); if(little) swabip(buf+412, 50*2);
00412 
00413   /* Write header */
00414   fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET);
00415   if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(4);
00416   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(5);
00417 
00418   return(0);
00419 }
00420 /*****************************************************************************/
00421 
00422 /*****************************************************************************/
00433 int ecat7WriteScanheader(FILE *fp, int blk, ECAT7_scanheader *h) {
00434   unsigned char buf[2*MatBLKSIZE];
00435   int little; /* 1 if current platform is little endian (i386), else 0 */
00436 
00437   if(ECAT7_TEST) printf("ecat7WriteScanheader()\n");
00438   if(fp==NULL || blk<2 || h==NULL) return(1);
00439   little=little_endian(); if(ECAT7_TEST) printf("little=%d\n", little);
00440   /* Clear buf */
00441   memset(buf, 0, 2*MatBLKSIZE);
00442   if(h->data_type==ECAT7_VAXI2) h->data_type=ECAT7_SUNI2;
00443   else if(h->data_type==ECAT7_VAXI4) h->data_type=ECAT7_SUNI4;
00444   else if(h->data_type==ECAT7_VAXR4) h->data_type=ECAT7_IEEER4;
00445 
00446   /* Copy the header fields and swap if necessary */
00447   memcpy(buf+0, &h->data_type, 2); if(little) swabip(buf+0, 2);
00448   memcpy(buf+2, &h->num_dimensions, 2); if(little) swabip(buf+2, 2);
00449   memcpy(buf+4, &h->num_r_elements, 2); if(little) swabip(buf+4, 2);
00450   memcpy(buf+6, &h->num_angles, 2); if(little) swabip(buf+6, 2);
00451   memcpy(buf+8, &h->corrections_applied, 2); if(little) swabip(buf+8, 2);
00452   memcpy(buf+10, h->num_z_elements, 64*2); if(little) swabip(buf+10, 64*2);
00453   memcpy(buf+138, &h->ring_difference, 2); if(little) swabip(buf+138, 2);
00454   memcpy(buf+140, &h->storage_order, 2); if(little) swabip(buf+140, 2);
00455   memcpy(buf+142, &h->axial_compression, 2); if(little) swabip(buf+142, 2);
00456   memcpy(buf+144, &h->x_resolution, 4); if(little) swawbip(buf+144, 4);
00457   memcpy(buf+148, &h->v_resolution, 4); if(little) swawbip(buf+148, 4);
00458   memcpy(buf+152, &h->z_resolution, 4); if(little) swawbip(buf+152, 4);
00459   memcpy(buf+156, &h->w_resolution, 4); if(little) swawbip(buf+156, 4);
00460   memcpy(buf+160, h->fill_gate, 6*2); if(little) swabip(buf+160, 6*2);
00461   memcpy(buf+172, &h->gate_duration, 4); if(little) swawbip(buf+172, 4);
00462   memcpy(buf+176, &h->r_wave_offset, 4); if(little) swawbip(buf+176, 4);
00463   memcpy(buf+180, &h->num_accepted_beats, 4); if(little) swawbip(buf+180, 4);
00464   memcpy(buf+184, &h->scale_factor, 4); if(little) swawbip(buf+184, 4);
00465   memcpy(buf+188, &h->scan_min, 2); if(little) swabip(buf+188, 2);
00466   memcpy(buf+190, &h->scan_max, 2); if(little) swabip(buf+190, 2);
00467   memcpy(buf+192, &h->prompts, 4); if(little) swawbip(buf+192, 4);
00468   memcpy(buf+196, &h->delayed, 4); if(little) swawbip(buf+196, 4);
00469   memcpy(buf+200, &h->multiples, 4); if(little) swawbip(buf+200, 4);
00470   memcpy(buf+204, &h->net_trues, 4); if(little) swawbip(buf+204, 4);
00471   memcpy(buf+208, &h->tot_avg_cor, 4); if(little) swawbip(buf+208, 4);
00472   memcpy(buf+212, &h->tot_avg_uncor, 4); if(little) swawbip(buf+212, 4);
00473   memcpy(buf+216, &h->total_coin_rate, 4); if(little) swawbip(buf+216, 4);
00474   memcpy(buf+220, &h->frame_start_time, 4); if(little) swawbip(buf+220, 4);
00475   memcpy(buf+224, &h->frame_duration, 4); if(little) swawbip(buf+224, 4);
00476   memcpy(buf+228, &h->deadtime_correction_factor, 4); if(little) swawbip(buf+228, 4);
00477   memcpy(buf+232, h->fill_cti, 90*2); if(little) swabip(buf+232, 90*2);
00478   memcpy(buf+412, h->fill_user, 50*2); if(little) swabip(buf+412, 50*2);
00479   memcpy(buf+512, h->uncor_singles, 128*4); if(little) swawbip(buf+512, 128*4);
00480 
00481   /* Write 3D scan header */
00482   fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(4);
00483   if(fwrite(buf, 1, 2*MatBLKSIZE, fp) != 2*MatBLKSIZE) return(5);
00484 
00485   return(0);
00486 }
00487 /*****************************************************************************/
00488 
00489 /*****************************************************************************/
00499 int ecat7Write2DScanheader(FILE *fp, int blk, ECAT7_2Dscanheader *h) {
00500   unsigned char buf[MatBLKSIZE];
00501   int little; /* 1 if current platform is little endian (i386), else 0 */
00502 
00503   if(ECAT7_TEST) printf("ecat7Write2DScanheader()\n");
00504   if(fp==NULL || blk<2 || h==NULL) return(1);
00505   little=little_endian(); if(ECAT7_TEST) printf("little=%d\n", little);
00506   /* Clear buf */
00507   memset(buf, 0, MatBLKSIZE);
00508   if(h->data_type==ECAT7_VAXI2) h->data_type=ECAT7_SUNI2;
00509   else if(h->data_type==ECAT7_VAXI4) h->data_type=ECAT7_SUNI4;
00510   else if(h->data_type==ECAT7_VAXR4) h->data_type=ECAT7_IEEER4;
00511 
00512   /* Copy the header fields and swap if necessary */
00513   memcpy(buf+0, &h->data_type, 2); if(little) swabip(buf+0, 2);
00514   memcpy(buf+2, &h->num_dimensions, 2); if(little) swabip(buf+2, 2);
00515   memcpy(buf+4, &h->num_r_elements, 2); if(little) swabip(buf+4, 2);
00516   memcpy(buf+6, &h->num_angles, 2); if(little) swabip(buf+6, 2);
00517   memcpy(buf+8, &h->corrections_applied, 2); if(little) swabip(buf+8, 2);
00518   memcpy(buf+10, &h->num_z_elements, 2); if(little) swabip(buf+10, 2);
00519   memcpy(buf+12, &h->ring_difference, 2); if(little) swabip(buf+12, 2);
00520   memcpy(buf+14, &h->x_resolution, 4); if(little) swawbip(buf+14, 4);
00521   memcpy(buf+18, &h->y_resolution, 4); if(little) swawbip(buf+18, 4);
00522   memcpy(buf+22, &h->z_resolution, 4); if(little) swawbip(buf+22, 4);
00523   memcpy(buf+26, &h->w_resolution, 4); if(little) swawbip(buf+26, 4);
00524   memcpy(buf+30, h->fill_gate, 6*2); if(little) swabip(buf+30, 6*2);
00525   memcpy(buf+42, &h->gate_duration, 4); if(little) swawbip(buf+42, 4);
00526   memcpy(buf+46, &h->r_wave_offset, 4); if(little) swawbip(buf+46, 4);
00527   memcpy(buf+50, &h->num_accepted_beats, 4); if(little) swawbip(buf+50, 4);
00528   memcpy(buf+54, &h->scale_factor, 4); if(little) swawbip(buf+54, 4);
00529   memcpy(buf+58, &h->scan_min, 2); if(little) swabip(buf+58, 2);
00530   memcpy(buf+60, &h->scan_max, 2); if(little) swabip(buf+60, 2);
00531   memcpy(buf+62, &h->prompts, 4); if(little) swawbip(buf+62, 4);
00532   memcpy(buf+66, &h->delayed, 4); if(little) swawbip(buf+66, 4);
00533   memcpy(buf+70, &h->multiples, 4); if(little) swawbip(buf+70, 4);
00534   memcpy(buf+74, &h->net_trues, 4); if(little) swawbip(buf+74, 4);
00535   memcpy(buf+78, h->cor_singles, 16*4); if(little) swawbip(buf+78, 16*4);
00536   memcpy(buf+142, h->uncor_singles, 16*4); if(little) swawbip(buf+142, 16*4);
00537   memcpy(buf+206, &h->tot_avg_cor, 4); if(little) swawbip(buf+206, 4);
00538   memcpy(buf+210, &h->tot_avg_uncor, 4); if(little) swawbip(buf+210, 4);
00539   memcpy(buf+214, &h->total_coin_rate, 4); if(little) swawbip(buf+214, 4);
00540   memcpy(buf+218, &h->frame_start_time, 4); if(little) swawbip(buf+218, 4);
00541   memcpy(buf+222, &h->frame_duration, 4); if(little) swawbip(buf+222, 4);
00542   memcpy(buf+226, &h->deadtime_correction_factor, 4); if(little) swawbip(buf+226, 4);
00543   memcpy(buf+230, h->physical_planes, 8*2); if(little) swabip(buf+230, 8*2);
00544   memcpy(buf+246, h->fill_cti, 83*2); if(little) swabip(buf+246, 83*2);
00545   memcpy(buf+412, h->fill_user, 50*2); if(little) swabip(buf+412, 50*2);
00546 
00547   /* Write header */
00548   fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET);
00549   if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(4);
00550   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(5);
00551 
00552   return(0);
00553 }
00554 /*****************************************************************************/
00555 
00556 /*****************************************************************************/
00566 int ecat7Write2DNormheader(FILE *fp, int blk, ECAT7_2Dnormheader *h) {
00567   unsigned char buf[MatBLKSIZE];
00568   int little; /* 1 if current platform is little endian (i386), else 0 */
00569 
00570   if(ECAT7_TEST) printf("ecat7Write2DNormheader()\n");
00571   if(fp==NULL || blk<2 || h==NULL) return(1);
00572   little=little_endian(); if(ECAT7_TEST) printf("little=%d\n", little);
00573   /* Clear buf */
00574   memset(buf, 0, MatBLKSIZE);
00575   if(h->data_type==ECAT7_VAXI2) h->data_type=ECAT7_SUNI2;
00576   else if(h->data_type==ECAT7_VAXI4) h->data_type=ECAT7_SUNI4;
00577   else if(h->data_type==ECAT7_VAXR4) h->data_type=ECAT7_IEEER4;
00578 
00579   /* Copy the header fields and swap if necessary */
00580   memcpy(buf+0, &h->data_type, 2); if(little) swabip(buf+0, 2);
00581   memcpy(buf+2, &h->num_dimensions, 2); if(little) swabip(buf+2, 2);
00582   memcpy(buf+4, &h->num_r_elements, 2); if(little) swabip(buf+4, 2);
00583   memcpy(buf+6, &h->num_angles, 2); if(little) swabip(buf+6, 2);
00584   memcpy(buf+8, &h->num_z_elements, 2); if(little) swabip(buf+8, 2);
00585   memcpy(buf+10, &h->ring_difference, 2); if(little) swabip(buf+10, 2);
00586   memcpy(buf+12, &h->scale_factor, 4); if(little) swawbip(buf+12, 4);
00587   memcpy(buf+16, &h->norm_min, 4); if(little) swawbip(buf+16, 4);
00588   memcpy(buf+20, &h->norm_max, 4); if(little) swawbip(buf+20, 4);
00589   memcpy(buf+24, &h->fov_source_width, 4); if(little) swawbip(buf+24, 4);
00590   memcpy(buf+28, &h->norm_quality_factor, 4); if(little) swawbip(buf+28, 4);
00591   memcpy(buf+32, &h->norm_quality_factor_code, 2); if(little) swabip(buf+32, 2);
00592   memcpy(buf+34, &h->storage_order, 2); if(little) swabip(buf+34, 2);
00593   memcpy(buf+36, &h->span, 2); if(little) swabip(buf+36, 2);
00594   memcpy(buf+38, h->fill_cti, 64*2); if(little) swabip(buf+38, 64*2);
00595   memcpy(buf+166, h->fill_cti, 123*2); if(little) swabip(buf+166, 123*2);
00596   memcpy(buf+412, h->fill_user, 50*2); if(little) swabip(buf+412, 50*2);
00597 
00598   /* Write header */
00599   fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET);
00600   if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(4);
00601   if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(5);
00602 
00603   return(0);
00604 }
00605 /*****************************************************************************/
00606 
00607 /*****************************************************************************/
00616 FILE *ecat7Create(const char *fname, ECAT7_mainheader *h) {
00617   FILE *fp;
00618   char tmp[FILENAME_MAX];
00619   int buf[MatBLKSIZE/4];
00620 
00621   if(ECAT7_TEST) printf("ecat7Create(%s, h)\n", fname);
00622   /* Check the arguments */
00623   if(fname==NULL || h==NULL) return(NULL);
00624   /* Check if file exists; backup, if necessary */
00625   if(access(fname, 0) != -1) {
00626     strcpy(tmp, fname); strcat(tmp, BACKUP_EXTENSION);
00627     if(access(tmp, 0) != -1) remove(tmp);
00628     if(ECAT7_TEST) printf("Renaming %s -> %s\n", fname, tmp);
00629     rename(fname, tmp);
00630   }
00631   /* Open file */
00632   fp=fopen(fname, "wb+"); if(fp==NULL) return(fp);
00633   /* Write main header */
00634   if(ecat7WriteMainheader(fp, h)) return(NULL);
00635   /* Construct an empty matrix list ; convert to little endian if necessary */
00636   memset(buf, 0, MatBLKSIZE);
00637   buf[0]=31; buf[1]=MatFirstDirBlk; if(little_endian()) swawbip(buf, MatBLKSIZE);
00638   /* Write data buffer */
00639   fseek(fp, (MatFirstDirBlk-1)*MatBLKSIZE, SEEK_SET);
00640   if(ftell(fp)!=(MatFirstDirBlk-1)*MatBLKSIZE) return(NULL);
00641   if(fwrite(buf, 4, MatBLKSIZE/4, fp) != MatBLKSIZE/4) return(NULL);
00642   /* OK, then return file pointer */
00643   return(fp);
00644 }
00645 /*****************************************************************************/
00646 
00647 /*****************************************************************************/
00657 int ecat7_is_scaling_needed(float amax, float *data, int nr) {
00658   int i;
00659   double d;
00660   
00661   if(nr<1 || data==NULL) return(0);
00662   /* scaling is necessary if all values are between -1 - 1 */
00663   if(amax<0.9999) return(1);
00664   /* Lets check first if at least the max value is close to integers or not */
00665   if(modf(amax, &d)>0.0001) return(1);
00666   /* if it is, then check all pixels */
00667   for(i=0; i<nr; i++) if(modf(*data++, &d)>0.0001) return(1);
00668   return(0);
00669 }
00670 /*****************************************************************************/
00671 
00672 /*****************************************************************************/
00682 int ecat7WriteImageMatrix(FILE *fp, int matrix_id, ECAT7_imageheader *h, float *fdata) {
00683   int i, nxtblk, blkNr, data_size, pxlNr, ret;
00684   float *fptr, fmin, fmax, g, f;
00685   char *mdata, *mptr;
00686   short int *sptr;
00687   
00688   
00689   if(ECAT7_TEST) printf("ecat7WriteImageMatrix(fp, %d, h, data)\n", matrix_id);
00690   if(fp==NULL || matrix_id<1 || h==NULL || fdata==NULL) {
00691     sprintf(ecat7errmsg, "invalid function parameter.\n");
00692     return(1);
00693   }
00694   if(h->data_type!=ECAT7_SUNI2) {
00695     sprintf(ecat7errmsg, "invalid data_type.\n");
00696     return(2);
00697   }
00698   /* nr of pixels */
00699   pxlNr=h->x_dimension*h->y_dimension;
00700   if(h->num_dimensions>2) pxlNr*=h->z_dimension;
00701   if(pxlNr<1) {
00702     sprintf(ecat7errmsg, "invalid matrix dimension.\n");
00703     return(3);
00704   }
00705   /* How much memory is needed for ALL pixels */
00706   data_size=pxlNr*ecat7pxlbytes(h->data_type);
00707   /* block nr taken by all pixels */
00708   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) {
00709     sprintf(ecat7errmsg, "invalid block number.\n");
00710     return(4);
00711   }
00712   /* Allocate memory for matrix data */
00713   mdata=(char*)calloc(blkNr, MatBLKSIZE); if(mdata==NULL) {
00714     sprintf(ecat7errmsg, "out of memory.\n");
00715     return(5);
00716   }
00717   /* Search for min and max for calculation of scale factor */
00718   fptr=fdata; fmin=fmax=*fptr;
00719   for(i=0; i<pxlNr; i++, fptr++) {
00720     if(*fptr>fmax) fmax=*fptr; else if(*fptr<fmin) fmin=*fptr;
00721   }
00722   if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
00723   if(g>0) f=32766./g; else f=1.0;
00724   /* Check if pixels values can be left as such with scale_factor = 1 */
00725   fptr=fdata;
00726   if(f>=1.0 && ecat7_is_scaling_needed(g, fptr, pxlNr)==0) f=1.0;
00727   /* Scale matrix data to shorts */
00728   h->scale_factor=1.0/f;
00729   sptr=(short int*)mdata; fptr=fdata;
00730   for(i=0; i<pxlNr; i++, sptr++, fptr++) *sptr=(short int)temp_roundf(f*(*fptr));
00731   /* Set header short min & max */
00732   h->image_min=(short int)temp_roundf(f*fmin);
00733   h->image_max=(short int)temp_roundf(f*fmax);
00734   /* Get block number for matrix header and data */
00735   nxtblk=ecat7EnterMatrix(fp, matrix_id, blkNr); if(nxtblk<1) {
00736     sprintf(ecat7errmsg, "cannot determine matrix block (%d).\n", -nxtblk);
00737     free(mdata); return(8);
00738   }
00739   if(ECAT7_TEST>2) printf("  block=%d fmin=%g fmax=%g\n", nxtblk, fmin, fmax);
00740   /* Write header */
00741   ret=ecat7WriteImageheader(fp, nxtblk, h); if(ret) {
00742     sprintf(ecat7errmsg, "cannot write subheader (%d).\n", ret);
00743     free(mdata); return(10);
00744   }
00745   /* Write matrix data */
00746   mptr=mdata;
00747   ret=ecat7WriteMatrixdata(fp, nxtblk+1, mptr, pxlNr, ecat7pxlbytes(h->data_type));
00748   free(mdata);
00749   if(ret) {
00750     sprintf(ecat7errmsg, "cannot write matrix data (%d).\n", ret);
00751     return(13);
00752   }
00753   return(0);
00754 }
00755 /*****************************************************************************/
00756 
00757 /*****************************************************************************/
00767 int ecat7Write2DScanMatrix(FILE *fp, int matrix_id, ECAT7_2Dscanheader *h, float *fdata) {
00768   int i, nxtblk, blkNr, data_size, pxlNr, ret;
00769   float *fptr, fmin, fmax, g, f;
00770   char *mdata, *mptr;
00771   short int *sptr;
00772 
00773 
00774   if(ECAT7_TEST) printf("ecat7Write2DScanMatrix(fp, %d, h, data)\n", matrix_id);
00775   if(fp==NULL || matrix_id<1 || h==NULL || fdata==NULL) {
00776     sprintf(ecat7errmsg, "invalid function parameter.\n");
00777     return(1);
00778   }
00779   if(h->data_type!=ECAT7_SUNI2) {
00780     sprintf(ecat7errmsg, "invalid data_type.\n");
00781     return(2);
00782   }
00783   /* nr of pixels */
00784   pxlNr=h->num_r_elements*h->num_angles;
00785   if(h->num_dimensions>2) pxlNr*=h->num_z_elements;
00786   if(pxlNr<1) {
00787     sprintf(ecat7errmsg, "invalid matrix dimension.\n");
00788     return(3);
00789   }
00790   /* How much memory is needed for ALL pixels */
00791   data_size=pxlNr*ecat7pxlbytes(h->data_type);
00792   /* block nr taken by all pixels */
00793   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) {
00794     sprintf(ecat7errmsg, "invalid block number.\n");
00795     return(4);
00796   }
00797   /* Allocate memory for matrix data */
00798   mdata=(char*)calloc(blkNr, MatBLKSIZE); if(mdata==NULL) {
00799     sprintf(ecat7errmsg, "out of memory.\n");
00800     return(5);
00801   }
00802   /* Search for min and max for calculation of scale factor */
00803   fptr=fdata; fmin=fmax=*fptr;
00804   for(i=0; i<pxlNr; i++, fptr++) {
00805     if(*fptr>fmax) fmax=*fptr; else if(*fptr<fmin) fmin=*fptr;
00806   }
00807   if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
00808   if(g>0) f=32766./g; else f=1.0;
00809   /* Check if pixels values can be left as such with scale_factor = 1 */
00810   fptr=fdata;
00811   if(f>=1.0 && ecat7_is_scaling_needed(g, fptr, pxlNr)==0) f=1.0;
00812   /* Scale matrix data to shorts */
00813   h->scale_factor=1.0/f;
00814   sptr=(short int*)mdata; fptr=fdata;
00815   for(i=0; i<pxlNr; i++, sptr++, fptr++) *sptr=(short int)temp_roundf(f*(*fptr));
00816   /* Set header short min & max */
00817   h->scan_min=(short int)temp_roundf(f*fmin);
00818   h->scan_max=(short int)temp_roundf(f*fmax);
00819   /* Get block number for matrix header and data */
00820   nxtblk=ecat7EnterMatrix(fp, matrix_id, blkNr); if(nxtblk<1) {
00821     sprintf(ecat7errmsg, "cannot determine matrix block (%d).\n", -nxtblk);
00822     free(mdata); return(8);
00823   }
00824   if(ECAT7_TEST>2) printf("  block=%d fmin=%g fmax=%g\n", nxtblk, fmin, fmax);
00825   /* Write header */
00826   ret=ecat7Write2DScanheader(fp, nxtblk, h); if(ret) {
00827     sprintf(ecat7errmsg, "cannot write subheader (%d).\n", ret);
00828     free(mdata); return(10);
00829   }
00830   /* Write matrix data */
00831   mptr=mdata;
00832   ret=ecat7WriteMatrixdata(fp, nxtblk+1, mptr, pxlNr, ecat7pxlbytes(h->data_type));
00833   free(mdata);
00834   if(ret) {
00835     sprintf(ecat7errmsg, "cannot write matrix data (%d).\n", ret);
00836     return(13);
00837   }
00838   return(0);
00839 }
00840 /*****************************************************************************/
00841 
00842 /*****************************************************************************/
00852 int ecat7WriteScanMatrix(FILE *fp, int matrix_id, ECAT7_scanheader *h, float *fdata) {
00853   int i, nxtblk, blkNr, data_size, pxlNr, dimz, ret;
00854   float *fptr, fmin, fmax, g, f;
00855   char *mdata, *mptr;
00856   short int *sptr;
00857 
00858 
00859   if(ECAT7_TEST) printf("ecat7WriteScanMatrix(fp, %d, h, data)\n", matrix_id);
00860   if(fp==NULL || matrix_id<1 || h==NULL || fdata==NULL) {
00861     sprintf(ecat7errmsg, "invalid function parameter.\n");
00862     return(1);
00863   }
00864   if(h->data_type!=ECAT7_SUNI2) {
00865     sprintf(ecat7errmsg, "invalid data_type.\n");
00866     return(2);
00867   }
00868   /* nr of pixels */
00869   pxlNr=h->num_r_elements*h->num_angles;
00870   for(i=dimz=0; i<64; i++) dimz+=h->num_z_elements[i]; pxlNr*=dimz;
00871   if(pxlNr<1) {
00872     sprintf(ecat7errmsg, "invalid matrix dimension.\n");
00873     return(3);
00874   }
00875   /* How much memory is needed for ALL pixels */
00876   data_size=pxlNr*ecat7pxlbytes(h->data_type);
00877   /* block nr taken by all pixels */
00878   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) {
00879     sprintf(ecat7errmsg, "invalid block number.\n");
00880     return(4);
00881   }
00882   /* Allocate memory for matrix data */
00883   mdata=(char*)calloc(blkNr, MatBLKSIZE); if(mdata==NULL) {
00884     sprintf(ecat7errmsg, "out of memory.\n");
00885     return(5);
00886   }
00887   /* Search for min and max for calculation of scale factor */
00888   fptr=fdata; fmin=fmax=*fptr;
00889   for(i=0; i<pxlNr; i++, fptr++) {
00890     if(*fptr>fmax) fmax=*fptr; else if(*fptr<fmin) fmin=*fptr;
00891   }
00892   if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
00893   if(g>0) f=32766./g; else f=1.0;
00894   /* Check if pixels values can be left as such with scale_factor = 1 */
00895   fptr=fdata;
00896   if(f>=1.0 && ecat7_is_scaling_needed(g, fptr, pxlNr)==0) f=1.0;
00897   /* Scale matrix data to shorts */
00898   h->scale_factor=1.0/f;
00899   sptr=(short int*)mdata; fptr=fdata;
00900   for(i=0; i<pxlNr; i++, sptr++, fptr++) *sptr=(short int)temp_roundf(f*(*fptr));
00901   /* Set header short min & max */
00902   h->scan_min=(short int)temp_roundf(f*fmin);
00903   h->scan_max=(short int)temp_roundf(f*fmax);
00904   /* Get block number for matrix header and data */
00905   /* Note that one extra block (blkNr+1) is needed for 3D scan header */
00906   nxtblk=ecat7EnterMatrix(fp, matrix_id, blkNr+1); if(nxtblk<1) {
00907     sprintf(ecat7errmsg, "cannot determine matrix block (%d).\n", -nxtblk);
00908     free(mdata); return(8);
00909   }
00910   if(ECAT7_TEST>2) printf("  block=%d fmin=%g fmax=%g\n", nxtblk, fmin, fmax);
00911   /* Write header */
00912   ret=ecat7WriteScanheader(fp, nxtblk, h); if(ret) {
00913     sprintf(ecat7errmsg, "cannot write subheader (%d).\n", ret);
00914     free(mdata); return(10);
00915   }
00916   /* Write matrix data */
00917   /* Note that 3D scan header takes TWO blocks */
00918   mptr=mdata;
00919   ret=ecat7WriteMatrixdata(fp, nxtblk+2, mptr, pxlNr, ecat7pxlbytes(h->data_type));
00920   free(mdata);
00921   if(ret) {
00922     sprintf(ecat7errmsg, "cannot write matrix data (%d).\n", ret);
00923     return(13);
00924   }
00925   return(0);
00926 }
00927 /*****************************************************************************/
00928 
00929 /*****************************************************************************/
00939 int ecat7WritePolarmapMatrix(FILE *fp, int matrix_id, ECAT7_polmapheader *h, float *fdata) {
00940   int i, nxtblk, blkNr, data_size, pxlNr, ret;
00941   float *fptr, fmin, fmax, g, f;
00942   char *mdata, *mptr;
00943   short int *sptr;
00944   
00945   
00946   if(ECAT7_TEST) printf("ecat7WritePolarmapMatrix(fp, %d, h, data)\n", matrix_id);
00947   if(fp==NULL || matrix_id<1 || h==NULL || fdata==NULL) {
00948     sprintf(ecat7errmsg, "invalid function parameter.\n");
00949     return(1);
00950   }
00951   if(h->data_type!=ECAT7_SUNI2) {
00952     sprintf(ecat7errmsg, "invalid data_type.\n");
00953     return(2);
00954   }
00955   /* nr of pixels */
00956   for(i=pxlNr=0; i<h->num_rings; i++) pxlNr+=h->sectors_per_ring[i];
00957   if(pxlNr<1) {
00958     sprintf(ecat7errmsg, "invalid matrix dimension.\n");
00959     return(3);
00960   }
00961   /* How much memory is needed for ALL pixels */
00962   data_size=pxlNr*ecat7pxlbytes(h->data_type);
00963   /* block nr taken by all pixels */
00964   blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) {
00965     sprintf(ecat7errmsg, "invalid block number.\n");
00966     return(4);
00967   }
00968   /* Allocate memory for matrix data */
00969   mdata=(char*)calloc(blkNr, MatBLKSIZE); if(mdata==NULL) {
00970     sprintf(ecat7errmsg, "out of memory.\n");
00971     return(5);
00972   }
00973   /* Search for min and max for calculation of scale factor */
00974   fptr=fdata; fmin=fmax=*fptr;
00975   for(i=0; i<pxlNr; i++, fptr++) {
00976     if(*fptr>fmax) fmax=*fptr; else if(*fptr<fmin) fmin=*fptr;
00977   }
00978   if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
00979   if(g>0) f=32766./g; else f=1.0;
00980   /* Check if pixels values can be left as such with scale_factor = 1 */
00981   fptr=fdata;
00982   if(f>=1.0 && ecat7_is_scaling_needed(g, fptr, pxlNr)==0) f=1.0;
00983   /* Scale matrix data to shorts */
00984   h->scale_factor=1.0/f;
00985   sptr=(short int*)mdata; fptr=fdata;
00986   for(i=0; i<pxlNr; i++, sptr++, fptr++) *sptr=(short int)temp_roundf(f*(*fptr));
00987   /* Set header short min & max */
00988   h->image_min=(short int)temp_roundf(f*fmin);
00989   h->image_max=(short int)temp_roundf(f*fmax);
00990   /* Get block number for matrix header and data */
00991   nxtblk=ecat7EnterMatrix(fp, matrix_id, blkNr); if(nxtblk<1) {
00992     sprintf(ecat7errmsg, "cannot determine matrix block (%d).\n", -nxtblk);
00993     free(mdata); return(8);
00994   }
00995   if(ECAT7_TEST>2) printf("  block=%d fmin=%g fmax=%g\n", nxtblk, fmin, fmax);
00996   /* Write header */
00997   ret=ecat7WritePolmapheader(fp, nxtblk, h); if(ret) {
00998     sprintf(ecat7errmsg, "cannot write subheader (%d).\n", ret);
00999     free(mdata); return(10);
01000   }
01001   /* Write matrix data */
01002   mptr=mdata;
01003   ret=ecat7WriteMatrixdata(fp, nxtblk+1, mptr, pxlNr, ecat7pxlbytes(h->data_type));
01004   free(mdata);
01005   if(ret) {
01006     sprintf(ecat7errmsg, "cannot write matrix data (%d).\n", ret);
01007     return(13);
01008   }
01009   return(0);
01010 }
01011 /*****************************************************************************/
01012 
01013 /*****************************************************************************/
01027 int ecat7WriteMatrixdata(FILE *fp, int start_block, char *data, int pxl_nr, int pxl_size) {
01028   unsigned char buf[MatBLKSIZE];
01029   char *dptr;
01030   int i, blkNr, dataSize, byteNr, little;
01031 
01032   if(ECAT7_TEST) printf("ecat7WriteMatrixdata(fp, %d, data, %d, %d)\n",
01033                         start_block, pxl_nr, pxl_size);
01034   if(fp==NULL || start_block<1 || data==NULL || pxl_nr<1 || pxl_size<1) return(1);
01035   little=little_endian(); memset(buf, 0, MatBLKSIZE);
01036   dataSize=pxl_nr*pxl_size;
01037   /* block nr taken by all pixels */
01038   blkNr=(dataSize+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(1);
01039   if(ECAT7_TEST>2) printf("    blkNr=%d\n", blkNr);
01040   /* Search the place for writing */
01041   fseek(fp, (start_block-1)*MatBLKSIZE, SEEK_SET);
01042   if(ftell(fp)!=(start_block-1)*MatBLKSIZE) return(2);
01043   /* Save blocks one at a time */
01044   for(i=0, dptr=data; i<blkNr && dataSize>0; i++) {
01045     byteNr=(dataSize<MatBLKSIZE)?dataSize:MatBLKSIZE;
01046     memcpy(buf, dptr, byteNr);
01047     /* Change matrix byte order in little endian platforms */
01048     if(little) {
01049       if(pxl_size==2) swabip(buf, byteNr);
01050       else if(pxl_size==4) swawbip(buf, byteNr);
01051     }
01052     /* Write block */
01053     if(fwrite(buf, 1, MatBLKSIZE, fp)!=MatBLKSIZE) return(3);
01054     /* Prepare for the next block */
01055     dptr+=byteNr; dataSize-=byteNr;
01056   } /* next block */
01057   return(0);
01058 }
01059 /*****************************************************************************/
01060 
01061 /*****************************************************************************/
01062