ecat63r.c
Go to the documentation of this file.
1 /******************************************************************************
2 
3  ecat63r.c (c) 2003,2004 by Turku PET Centre
4 
5  Procedures for reading ECAT 6.3 format.
6 
7  Assumptions:
8  1. Data_type in headers specifies, whether ints, long ints and floats in
9  header and matrix data are in VAX format (1, 2, 3 and 4) or in IEEE
10  (i386, SUN) format.
11  2. Data is automatically converted to little or big endian when read,
12  according to the current platform.
13  3. Data is automatically converted out from the VAX format when read,
14  but header data_type remains as original.
15 
16  Version: Contents were previously in ecat63.c
17  2003-07-21 Vesa Oikonen
18  2003-09-05 VO
19  Added function ecat63ReadImageMatrix() and ecat63ReadScanMatrix().
20  2003-09-08 VO
21  Added function ecat63pxlbytes().
22  2004-06-21 VO
23  Larger range allowed for calibration factor in ecat63ReadMainheader().
24  2004-09-20 VO
25  Doxygen style comments.
26 
27 ******************************************************************************/
28 #include <stdio.h>
29 #include <stdlib.h>
30 #include <math.h>
31 #include <ctype.h>
32 #include <string.h>
33 #include <unistd.h>
34 #include <time.h>
35 /*****************************************************************************/
36 #include "swap.h"
37 #include "include/ecat63.h"
38 /*****************************************************************************/
39 
40 /*****************************************************************************/
51  unsigned char buf[MatBLKSIZE];
52  int i;
53  int little; /* 1 if current platform is little endian (i386), else 0 */
54  int vaxdata=1; /* 1 if data is in VAX format, else 0 */
55 
56  if(ECAT63_TEST) printf("ecat63ReadMainheader()\n");
57  if(fp==NULL || h==NULL) return(1);
58  little=little_endian();
59  /* Seek the first block */
60  fseek(fp, 0, SEEK_SET); if(ftell(fp)!=0) return(1);
61  /* Read the subheader block */
62  if(fread(buf, MatBLKSIZE, 1, fp)<1) return(2);
63 
64  /* Copy char data to header structure */
65  memcpy(h->ecat_format, buf+0, 14); memcpy(h->fill1, buf+14, 14);
66  memcpy(h->original_file_name, buf+28, 20); memcpy(h->node_id, buf+56, 10);
67  memcpy(h->isotope_code, buf+78, 8); memcpy(h->radiopharmaceutical, buf+90, 32);
68  memcpy(h->study_name, buf+162, 12); memcpy(h->patient_id, buf+174, 16);
69  memcpy(h->patient_name, buf+190, 32); h->patient_sex=buf[222];
70  memcpy(h->patient_age, buf+223, 10); memcpy(h->patient_height, buf+233, 10);
71  memcpy(h->patient_weight, buf+243, 10); h->patient_dexterity=buf[253];
72  memcpy(h->physician_name, buf+254, 32); memcpy(h->operator_name, buf+286, 32);
73  memcpy(h->study_description, buf+318, 32); memcpy(h->facility_name, buf+356, 20);
74  memcpy(h->user_process_code, buf+462, 10);
75 
76  /* Copy short ints */
77  /* in big endian platform, change byte order temporarily */
78  if(!little) swabip(buf, MatBLKSIZE);
79  memcpy(&h->data_type, buf+50, 2); if(h->data_type>4) vaxdata=0;
80  /*printf("main_header.data_type=%d\n", h->data_type);*/
81  memcpy(&h->sw_version, buf+48, 2);
82  memcpy(&h->system_type, buf+52, 2); memcpy(&h->file_type, buf+54, 2);
83  memcpy(&h->scan_start_day, buf+66, 2); memcpy(&h->scan_start_month, buf+68, 2);
84  memcpy(&h->scan_start_year, buf+70, 2); memcpy(&h->scan_start_hour, buf+72, 2);
85  memcpy(&h->scan_start_minute, buf+74, 2); memcpy(&h->scan_start_second, buf+76, 2);
86  memcpy(&h->rot_source_speed, buf+134, 2); memcpy(&h->wobble_speed, buf+136, 2);
87  memcpy(&h->transm_source_type, buf+138, 2); memcpy(&h->transaxial_samp_mode, buf+148, 2);
88  memcpy(&h->coin_samp_mode, buf+150, 2); memcpy(&h->axial_samp_mode, buf+152, 2);
89  memcpy(&h->calibration_units, buf+158, 2); memcpy(&h->compression_code, buf+160, 2);
90  memcpy(&h->acquisition_type, buf+350, 2); memcpy(&h->bed_type, buf+352, 2);
91  memcpy(&h->septa_type, buf+354, 2); memcpy(&h->num_planes, buf+376, 2);
92  memcpy(&h->num_frames, buf+378, 2); memcpy(&h->num_gates, buf+380, 2);
93  memcpy(&h->num_bed_pos, buf+382, 2); memcpy(&h->lwr_sctr_thres, buf+452, 2);
94  memcpy(&h->lwr_true_thres, buf+454, 2); memcpy(&h->upr_true_thres, buf+456, 2);
95  memcpy(h->fill2, buf+472, 40);
96  /* Change back the byte order */
97  if(!little) swabip(buf, MatBLKSIZE);
98 
99  /* Copy floats */
100  h->isotope_halflife=ecat63rFloat(buf+86, vaxdata, little);
101  h->gantry_tilt=ecat63rFloat(buf+122, vaxdata, little);
102  h->gantry_rotation=ecat63rFloat(buf+126, vaxdata, little);
103  h->bed_elevation=ecat63rFloat(buf+130, vaxdata, little);
104  h->axial_fov=ecat63rFloat(buf+140, vaxdata, little);
105  h->transaxial_fov=ecat63rFloat(buf+144, vaxdata, little);
106  h->calibration_factor=ecat63rFloat(buf+154, vaxdata, little);
107  h->init_bed_position=ecat63rFloat(buf+384, vaxdata, little);
108  for(i=0; i<15; i++) h->bed_offset[i]=ecat63rFloat(buf+388+i*4, vaxdata, little);
109  h->plane_separation=ecat63rFloat(buf+448, vaxdata, little);
110  h->collimator=ecat63rFloat(buf+458, vaxdata, little);
111 
112  /* Check file format and platform */
113  if(ECAT63_TEST) printf("ecat_format='%.14s'\n", h->ecat_format);
114  /* if format is not specified, ECAT63 is assumed */
115  if(h->ecat_format[0]==(char)0) {
116  strcpy(h->ecat_format, "ECAT63");
117  }
118  /* only ECAT63 format is approved here */
119  if(ECAT63_TEST) printf("ecat_format='%.14s'\n", h->ecat_format);
120  if(strncmp(h->ecat_format, "ECAT63", 6)!=0) return(3);
121 
122  /* Check that most important contents are ok */
123  if(h->data_type<BYTE_TYPE || h->data_type>SUN_I4) {
124  if(ECAT63_TEST) printf("Invalid data types; probable conversion error.\n");
125  return(5);
126  }
127  if(h->calibration_factor<0.0 || h->calibration_factor>1.0e12) {
128  if(ECAT63_TEST) printf("Invalid calibration factor; possible conversion error.\n");
129  return(6);
130  }
131  if(h->file_type!=RAW_DATA && h->file_type!=IMAGE_DATA &&
132  h->file_type!=ATTN_DATA && h->file_type!=NORM_DATA) {
133  if(ECAT63_TEST) printf("Invalid file types; probable conversion error.\n");
134  return(7);
135  }
136 
137  return(0);
138 }
139 /*****************************************************************************/
140 
141 /*****************************************************************************/
152 int ecat63ReadImageheader(FILE *fp, int blk, ECAT63_imageheader *h) {
153  unsigned char buf[MatBLKSIZE];
154  int i;
155  int little; /* 1 if current platform is little endian (i386), else 0 */
156  int vaxdata=1; /* 1 if data is in VAX format, else 0 */
157 
158  if(ECAT63_TEST) printf("ecat63ReadImageheader(fp, %d ih)\n", blk);
159  if(fp==NULL || blk<2 || h==NULL) return(1);
160  little=little_endian();
161  /* Seek the subheader block */
162  fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(2);
163  /* Read the subheader block */
164  if(fread(buf, MatBLKSIZE, 1, fp)<1) return(3);
165 
166  /* Copy char data to header structure */
167  memcpy(h->fill1, buf+0, 126); memcpy(h->annotation, buf+420, 40);
168 
169  /* Copy short ints */
170  /* in big endian platform, change byte order temporarily */
171  if(!little) swabip(buf, MatBLKSIZE);
172  memcpy(&h->data_type, buf+126, 2); if(h->data_type>4) vaxdata=0;
173  /*printf("data_type=%d\n", h->data_type);*/
174  memcpy(&h->num_dimensions, buf+128, 2);
175  memcpy(&h->dimension_1, buf+132, 2); memcpy(&h->dimension_2, buf+134, 2);
176  memcpy(&h->image_min, buf+176, 2); memcpy(&h->image_max, buf+178, 2);
177  memcpy(&h->slice_location, buf+200, 2); memcpy(&h->recon_start_hour, buf+202, 2);
178  memcpy(&h->recon_start_min, buf+204, 2); memcpy(&h->recon_start_sec, buf+206, 2);
179  memcpy(&h->filter_code, buf+236, 2); memcpy(&h->processing_code, buf+376, 2);
180  memcpy(&h->quant_units, buf+380, 2); memcpy(&h->recon_start_day, buf+382, 2);
181  memcpy(&h->recon_start_month, buf+384, 2); memcpy(&h->recon_start_year, buf+386, 2);
182  memcpy(h->fill2, buf+460, 52);
183  /* Change back the byte order */
184  if(!little) swabip(buf, MatBLKSIZE);
185 
186  /* Copy ints */
187  h->frame_duration=ecat63rInt(buf+192, vaxdata, little);
188  h->frame_start_time=ecat63rInt(buf+196, vaxdata, little);
189  h->recon_duration=ecat63rInt(buf+208, vaxdata, little);
190  h->scan_matrix_num=ecat63rInt(buf+238, vaxdata, little);
191  h->norm_matrix_num=ecat63rInt(buf+242, vaxdata, little);
192  h->atten_cor_mat_num=ecat63rInt(buf+246, vaxdata, little);
193 
194  /* Copy floats */
195  h->x_origin=ecat63rFloat(buf+160, vaxdata, little);
196  h->y_origin=ecat63rFloat(buf+164, vaxdata, little);
197  h->recon_scale=ecat63rFloat(buf+168, vaxdata, little);
198  h->quant_scale=ecat63rFloat(buf+172, vaxdata, little);
199  h->pixel_size=ecat63rFloat(buf+184, vaxdata, little);
200  h->slice_width=ecat63rFloat(buf+188, vaxdata, little);
201  h->image_rotation=ecat63rFloat(buf+296, vaxdata, little);
202  h->plane_eff_corr_fctr=ecat63rFloat(buf+300, vaxdata, little);
203  h->decay_corr_fctr=ecat63rFloat(buf+304, vaxdata, little);
204  h->loss_corr_fctr=ecat63rFloat(buf+308, vaxdata, little);
205  h->intrinsic_tilt=ecat63rFloat(buf+312, vaxdata, little);
206  h->ecat_calibration_fctr=ecat63rFloat(buf+388, vaxdata, little);
207  h->well_counter_cal_fctr=ecat63rFloat(buf+392, vaxdata, little);
208  for(i=0; i<6; i++) h->filter_params[i]=ecat63rFloat(buf+396+i*4, vaxdata, little);
209 
210  /* Check that most important contents are ok */
211  if(h->data_type<BYTE_TYPE || h->data_type>SUN_I4) {
212  if(ECAT63_TEST) printf("Invalid data types; probable conversion error.\n");
213  return(4);
214  }
215  if(h->ecat_calibration_fctr<0.0 || h->ecat_calibration_fctr>1.0e10) {
216  if(ECAT63_TEST) printf("Invalid calibration factor; probable conversion error.\n");
217  return(5);
218  }
219  if(h->frame_duration<0.0 || h->frame_duration>1.0e12) {
220  if(ECAT63_TEST) printf("Invalid frame duration; probable conversion error.\n");
221  return(6);
222  }
223 
224  return(0);
225 }
226 /*****************************************************************************/
227 
228 /*****************************************************************************/
238 int ecat63ReadAttnheader(FILE *fp, int blk, ECAT63_attnheader *h) {
239  unsigned char buf[MatBLKSIZE];
240  int little; /* 1 if current platform is little endian (i386), else 0 */
241  int vaxdata=1; /* 1 if data is in VAX format, else 0 */
242 
243  if(ECAT63_TEST) printf("ecat63ReadAttnheader(fp, %d, ah)\n", blk);
244  if(fp==NULL || blk<2 || h==NULL) return(1);
245  little=little_endian();
246  /* Seek the subheader block */
247  fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(2);
248  /* Read the subheader block */
249  if(fread(buf, MatBLKSIZE, 1, fp)<1) return(3);
250 
251  /* Copy short ints */
252  /* in big endian platform, change byte order temporarily */
253  if(!little) swabip(buf, MatBLKSIZE);
254  memcpy(&h->data_type, buf+126, 2); if(h->data_type>4) vaxdata=0;
255  /*printf("data_type=%d\n", h->data_type);*/
256  memcpy(&h->attenuation_type, buf+128, 2);
257  memcpy(&h->dimension_1, buf+132, 2); memcpy(&h->dimension_2, buf+134, 2);
258  /* Change back the byte order */
259  if(!little) swabip(buf, MatBLKSIZE);
260 
261  /* Copy floats */
262  h->scale_factor=ecat63rFloat(buf+182, vaxdata, little);
263  h->x_origin=ecat63rFloat(buf+186, vaxdata, little);
264  h->y_origin=ecat63rFloat(buf+190, vaxdata, little);
265  h->x_radius=ecat63rFloat(buf+194, vaxdata, little);
266  h->y_radius=ecat63rFloat(buf+198, vaxdata, little);
267  h->tilt_angle=ecat63rFloat(buf+202, vaxdata, little);
268  h->attenuation_coeff=ecat63rFloat(buf+206, vaxdata, little);
269  h->sample_distance=ecat63rFloat(buf+210, vaxdata, little);
270 
271  /* Check that most important contents are ok */
272  if(h->data_type<BYTE_TYPE || h->data_type>SUN_I4) {
273  if(ECAT63_TEST) printf("Invalid data types; probable conversion error.\n");
274  return(4);
275  }
276  if(h->scale_factor<=0.0 || h->scale_factor>1.0e8) {
277  if(ECAT63_TEST) printf("Invalid scale factor; probable conversion error.\n");
278  return(5);
279  }
280 
281  return(0);
282 }
283 /*****************************************************************************/
284 
285 /*****************************************************************************/
296 int ecat63ReadScanheader(FILE *fp, int blk, ECAT63_scanheader *h) {
297  unsigned char buf[MatBLKSIZE];
298  int i;
299  int little; /* 1 if current platform is little endian (i386), else 0 */
300  int vaxdata=1; /* 1 if data is in VAX format, else 0 */
301 
302  if(ECAT63_TEST) printf("ecat63ReadScanheader(fp, %d, sh)\n", blk);
303  if(fp==NULL || blk<2 || h==NULL) return(1);
304  little=little_endian();
305  /* Seek the subheader block */
306  fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(2);
307  /* Read the subheader block */
308  if(fread(buf, MatBLKSIZE, 1, fp)<1) return(3);
309 
310  /* Copy char data to header structure */
311  memcpy(h->fill1, buf+0, 126);
312 
313  /* Copy short ints */
314  /* in big endian platform, change byte order temporarily */
315  if(!little) swabip(buf, MatBLKSIZE);
316  memcpy(&h->data_type, buf+126, 2); if(h->data_type>4) vaxdata=0;
317  /*printf("data_type=%d\n", h->data_type);*/
318  memcpy(&h->dimension_1, buf+132, 2); memcpy(&h->dimension_2, buf+134, 2);
319  memcpy(&h->smoothing, buf+136, 2); memcpy(&h->processing_code, buf+138, 2);
320  memcpy(&h->frame_duration_sec, buf+170, 2);
321  memcpy(&h->scan_min, buf+192, 2); memcpy(&h->scan_max, buf+194, 2);
322  memcpy(h->fill2, buf+468, 44);
323  /* Change back the byte order */
324  if(!little) swabip(buf, MatBLKSIZE);
325 
326  /* Copy ints */
327  h->gate_duration=ecat63rInt(buf+172, vaxdata, little);
328  h->r_wave_offset=ecat63rInt(buf+176, vaxdata, little);
329  h->prompts=ecat63rInt(buf+196, vaxdata, little);
330  h->delayed=ecat63rInt(buf+200, vaxdata, little);
331  h->multiples=ecat63rInt(buf+204, vaxdata, little);
332  h->net_trues=ecat63rInt(buf+208, vaxdata, little);
333  h->total_coin_rate=ecat63rInt(buf+452, vaxdata, little);
334  h->frame_start_time=ecat63rInt(buf+456, vaxdata, little);
335  h->frame_duration=ecat63rInt(buf+460, vaxdata, little);
336 
337  /* Copy floats */
338  h->sample_distance=ecat63rFloat(buf+146, vaxdata, little);
339  h->isotope_halflife=ecat63rFloat(buf+166, vaxdata, little);
340  h->scale_factor=ecat63rFloat(buf+182, vaxdata, little);
341  for(i=0; i<16; i++) h->cor_singles[i]=ecat63rFloat(buf+316+i*4, vaxdata, little);
342  for(i=0; i<16; i++) h->uncor_singles[i]=ecat63rFloat(buf+380+i*4, vaxdata, little);
343  h->tot_avg_cor=ecat63rFloat(buf+444, vaxdata, little);
344  h->tot_avg_uncor=ecat63rFloat(buf+448, vaxdata, little);
345  h->loss_correction_fctr=ecat63rFloat(buf+464, vaxdata, little);
346 
347  /* Check that most important contents are ok */
348  if(h->data_type<BYTE_TYPE || h->data_type>SUN_I4) {
349  if(ECAT63_TEST) printf("Invalid data types; probable conversion error.\n");
350  return(4);
351  }
352  if(h->scale_factor<=0.0 || h->scale_factor>1.0e8) {
353  if(ECAT63_TEST) printf("Invalid scale factor; probable conversion error.\n");
354  return(5);
355  }
356  if(h->frame_duration<0.0 || h->frame_duration>1.0e12) {
357  if(ECAT63_TEST) printf("Invalid frame duration; probable conversion error.\n");
358  return(6);
359  }
360 
361  return(0);
362 }
363 /*****************************************************************************/
364 
365 /*****************************************************************************/
375 int ecat63ReadNormheader(FILE *fp, int blk, ECAT63_normheader *h) {
376  unsigned char buf[MatBLKSIZE];
377  int little; /* 1 if current platform is little endian (i386), else 0 */
378  int vaxdata=1; /* 1 if data is in VAX format, else 0 */
379 
380  if(ECAT63_TEST) printf("ecat63ReadNormheader(fp, %d, nh)\n", blk);
381  if(fp==NULL || blk<2 || h==NULL) return(1);
382  little=little_endian();
383  /* Seek the subheader block */
384  fseek(fp, (blk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(blk-1)*MatBLKSIZE) return(2);
385  /* Read the subheader block */
386  if(fread(buf, MatBLKSIZE, 1, fp)<1) return(3);
387 
388  /* Copy short ints */
389  /* in big endian platform, change byte order temporarily */
390  if(!little) swabip(buf, MatBLKSIZE);
391  memcpy(&h->data_type, buf+126, 2); if(h->data_type>4) vaxdata=0;
392  /*printf("data_type=%d\n", h->data_type);*/
393  memcpy(&h->dimension_1, buf+132, 2); memcpy(&h->dimension_2, buf+134, 2);
394  memcpy(&h->norm_hour, buf+186, 2); memcpy(&h->norm_minute, buf+188, 2);
395  memcpy(&h->norm_second, buf+190, 2); memcpy(&h->norm_day, buf+192, 2);
396  memcpy(&h->norm_month, buf+194, 2); memcpy(&h->norm_year, buf+196, 2);
397  /* Change back the byte order */
398  if(!little) swabip(buf, MatBLKSIZE);
399 
400  /* Copy floats */
401  h->scale_factor=ecat63rFloat(buf+182, vaxdata, little);
402  h->fov_source_width=ecat63rFloat(buf+198, vaxdata, little);
403 
404  /* Check that most important contents are ok */
405  if(h->data_type<BYTE_TYPE && h->data_type>SUN_I4) {
406  if(ECAT63_TEST) printf("Invalid data types; probable conversion error.\n");
407  return(4);
408  }
409  if(h->scale_factor<=0.0 || h->scale_factor>1.0e8) {
410  if(ECAT63_TEST) printf("Invalid scale factor; probable conversion error.\n");
411  return(5);
412  }
413 
414  return(0);
415 }
416 /*****************************************************************************/
417 
418 /*****************************************************************************/
432 int ecat63ReadMatdata(FILE *fp, int strtblk, int blkNr, char *data, int dtype) {
433  int i, n, little, err=0;
434  char *cptr;
435  float f;
436 
437 
438  if(ECAT63_TEST) printf("ecat63ReadMatdata(fp, %d, %d, data, %d)\n", strtblk, blkNr, dtype);
439  /* Check the arguments */
440  if(blkNr<=0 || strtblk<1 || data==NULL) return(1);
441  /* Seek the first data block */
442  fseek(fp, (strtblk-1)*MatBLKSIZE, SEEK_SET);
443  if(ftell(fp)!=(strtblk-1)*MatBLKSIZE) return(9);
444  /* Read the data blocks */
445  if(fread(data, MatBLKSIZE, blkNr, fp) < blkNr) return(2);
446  /* Translate data if necessary */
447  little=little_endian();
448  switch(dtype) {
449  case BYTE_TYPE: /* byte format...no translation necessary */
450  break;
451  case VAX_I2: /* byte conversion necessary on big endian platform */
452  if(!little) {cptr=data; swabip(cptr, blkNr*MatBLKSIZE);}
453  break;
454  case VAX_I4:
455  for(i=0, cptr=data; i<blkNr*MatBLKSIZE; i+=4, cptr+=4) {
456  n=ecat63rInt(cptr, 1, little); memcpy(cptr, &n, 4);
457  }
458  break;
459  case VAX_R4:
460  for(i=0, cptr=data; i<blkNr*MatBLKSIZE; i+=4, cptr+=4) {
461  f=ecat63rFloat(cptr, 1, little); memcpy(cptr, &f, 4);
462  }
463  break;
464  case IEEE_R4: /* IEEE float ; byte conversion necessary on big end platforms */
465  case SUN_I4: /* SUN int ; byte conversion necessary on big end platforms */
466  if(!little) swawbip(data, blkNr*MatBLKSIZE);
467  break;
468  case SUN_I2: /* SUN short ; byte conversion necessary on big end platforms */
469  if(!little) swabip(data, blkNr*MatBLKSIZE);
470  break;
471  default: /* if something else, for now think it as an error */
472  err=2;
473  break;
474  }
475  return(err);
476 }
477 /*****************************************************************************/
478 
479 /*****************************************************************************/
494 int ecat63ReadImageMatrix(FILE *fp, int first_block, int last_block, ECAT63_imageheader *h, float **fdata) {
495  int i, ret, blockNr, pxlNr;
496  char *mdata, *mptr;
497  float *_fdata, *fptr;
498  short int *sptr;
499  int *iptr;
500 
501 
502  if(ECAT63_TEST) printf("ecat63ReadImageMatrix(fp, %d, %d, hdr, fdata)\n",
503  first_block, last_block);
504  if(fp==NULL || first_block<=MatFirstDirBlk || h==NULL) {
505  sprintf(ecat63errmsg, "invalid function parameter.\n");
506  return(1);
507  }
508  *fdata=(float*)NULL;
509 
510  /* Read subheader */
511  ret=ecat63ReadImageheader(fp, first_block, h);
512  if(ret) {
513  sprintf(ecat63errmsg, "cannot read subheader (%d).\n", ret);
514  return(5);
515  }
516  if(ECAT63_TEST>4) ecat63PrintImageheader(h, stdout);
517  pxlNr=h->dimension_1*h->dimension_2;
518  if(pxlNr<=0) {
519  sprintf(ecat63errmsg, "invalid matrix dimension.\n");
520  return(6);
521  }
522 
523  /* Read matrix data */
524  blockNr=last_block-first_block; if(blockNr<1) return(0);
525  mdata=(char*)malloc(blockNr*MatBLKSIZE);
526  if(mdata==NULL) {
527  sprintf(ecat63errmsg, "cannot allocate memory.\n");
528  return(8);
529  }
530  mptr=mdata;
531  ret=ecat63ReadMatdata(fp, first_block+1, blockNr, mptr, h->data_type);
532  if(ret || mdata==NULL) {
533  sprintf(ecat63errmsg, "cannot read matrix data (%d).\n", ret);
534  free(mdata); return(9);
535  }
536 
537  /* Allocate memory for float data */
538  _fdata=(float*)malloc(pxlNr*sizeof(float));
539  if(_fdata==NULL) {
540  sprintf(ecat63errmsg, "cannot allocate memory.\n");
541  free(mdata); return(11);
542  }
543 
544  /* Convert matrix data to floats */
546  fptr=_fdata; mptr=mdata;
547  if(h->data_type==BYTE_TYPE) {
548  for(i=0; i<pxlNr; i++, mptr++, fptr++)
549  *fptr=h->quant_scale*(float)(*mptr);
550  } else if(h->data_type==VAX_I2 || h->data_type==SUN_I2) {
551  for(i=0; i<pxlNr; i++, mptr+=2, fptr++) {
552  sptr=(short int*)mptr;
553  *fptr=h->quant_scale*(float)(*sptr);
554  }
555  } else if(h->data_type==VAX_I4 || h->data_type==SUN_I4) {
556  for(i=0; i<pxlNr; i++, mptr+=4, fptr++) {
557  iptr=(int*)mptr;
558  *fptr=h->quant_scale*(float)(*iptr);
559  }
560  } else if(h->data_type==VAX_R4 || h->data_type==IEEE_R4) {
561  memcpy(fptr, mptr, pxlNr*4);
562  for(i=0; i<pxlNr; i++, fptr++) *fptr *= h->quant_scale;
563  }
564  free(mdata);
565  *fdata=_fdata;
566 
567  return(0);
568 }
569 /*****************************************************************************/
570 
571 /*****************************************************************************/
586 int ecat63ReadScanMatrix(FILE *fp, int first_block, int last_block, ECAT63_scanheader *h, float **fdata) {
587  int i, ret, blockNr, pxlNr;
588  char *mdata, *mptr;
589  float *_fdata, *fptr;
590  short int *sptr;
591  int *iptr;
592 
593 
594  if(ECAT63_TEST) printf("ecat63ReadScanMatrix(fp, %d, %d, hdr, fdata)\n",
595  first_block, last_block);
596  if(fp==NULL || first_block<=MatFirstDirBlk || h==NULL) {
597  sprintf(ecat63errmsg, "invalid function parameter.\n");
598  return(1);
599  }
600  *fdata=(float*)NULL;
601 
602  /* Read subheader */
603  ret=ecat63ReadScanheader(fp, first_block, h);
604  if(ret) {
605  sprintf(ecat63errmsg, "cannot read subheader (%d).\n", ret);
606  return(5);
607  }
608  if(ECAT63_TEST>4) ecat63PrintScanheader(h, stdout);
609  pxlNr=h->dimension_1*h->dimension_2;
610  if(pxlNr<=0) {
611  sprintf(ecat63errmsg, "invalid matrix dimension.\n");
612  return(6);
613  }
614 
615  /* Read matrix data */
616  blockNr=last_block-first_block; if(blockNr<1) return(0);
617  mdata=(char*)malloc(blockNr*MatBLKSIZE);
618  if(mdata==NULL) {
619  sprintf(ecat63errmsg, "cannot allocate memory.\n");
620  return(8);
621  }
622  mptr=mdata;
623  ret=ecat63ReadMatdata(fp, first_block+1, blockNr, mptr, h->data_type);
624  if(ret || mdata==NULL) {
625  sprintf(ecat63errmsg, "cannot read matrix data (%d).\n", ret);
626  free(mdata); return(9);
627  }
628 
629  /* Allocate memory for float data */
630  _fdata=(float*)malloc(pxlNr*sizeof(float));
631  if(_fdata==NULL) {
632  sprintf(ecat63errmsg, "cannot allocate memory.\n");
633  free(mdata); return(11);
634  }
635 
636  /* Convert matrix data to floats */
637  fptr=_fdata; mptr=mdata;
638  if(h->data_type==BYTE_TYPE) {
639  for(i=0; i<pxlNr; i++, mptr++, fptr++)
640  *fptr=h->scale_factor*(float)(*mptr);
641  } else if(h->data_type==VAX_I2 || h->data_type==SUN_I2) {
642  for(i=0; i<pxlNr; i++, mptr+=2, fptr++) {
643  sptr=(short int*)mptr;
644  *fptr=h->scale_factor*(float)(*sptr);
645  }
646  } else if(h->data_type==VAX_I4 || h->data_type==SUN_I4) {
647  for(i=0; i<pxlNr; i++, mptr+=4, fptr++) {
648  iptr=(int*)mptr;
649  *fptr=h->scale_factor*(float)(*iptr);
650  }
651  } else if(h->data_type==VAX_R4 || h->data_type==IEEE_R4) {
652  memcpy(fptr, mptr, pxlNr*4);
653  for(i=0; i<pxlNr; i++, fptr++) *fptr *= h->scale_factor;
654  }
655  free(mdata);
656  *fdata=_fdata;
657 
658  return(0);
659 }
660 /*****************************************************************************/
661 
662 /*****************************************************************************/
671 float ecat63rFloat(void *bufi, int isvax, int islittle) {
672  union {unsigned int ul; float f;} t;
673 
674  memcpy(&t.ul, bufi, 4); if(t.ul==0) {return(0.0);}
675  if(isvax) { /* if input is in VAX format */
676  /* Swap words on i386 and bytes on SUN */
677  if(islittle) swawip(&t.ul, 4); else swabip(&t.ul, 4);
678  t.ul-=(2L<<23); /* subtract 2 from exp */
679  } else { /* input is in i386 format */
680  if(!islittle) swawbip(&t.ul, 4); /* Switch words and bytes on SUN */
681  }
682  return(t.f);
683 }
684 
694 int ecat63rInt(void *bufi, int isvax, int islittle) {
695  int i;
696 
697  /* Swap both words and bytes on SUN */
698  memcpy(&i, bufi, 4); if(!islittle) swawbip(&i, 4);
699  return(i);
700 }
701 /*****************************************************************************/
702 
703 /*****************************************************************************/
711 int ecat63pxlbytes(short int data_type) {
712  int byteNr=0;
713  switch(data_type) {
714  case BYTE_TYPE: byteNr=1; break;
715  case VAX_I2:
716  case SUN_I2: byteNr=2; break;
717  case VAX_I4:
718  case VAX_R4:
719  case IEEE_R4:
720  case SUN_I4: byteNr=4; break;
721  }
722  return(byteNr);
723 }
724 /*****************************************************************************/
725 
726 /*****************************************************************************/
727