ecat63w.c
Go to the documentation of this file.
1 /******************************************************************************
2 
3  ecat63w.c (c) 2003-2008 by Turku PET Centre
4 
5  Procedures for writing ECAT 6.3 matrix data.
6 
7  This program is free software; you can redistribute it and/or modify it under
8  the terms of the GNU General Public License as published by the Free Software
9  Foundation; either version 2 of the License, or (at your option) any later
10  version.
11 
12  This program is distributed in the hope that it will be useful, but WITHOUT
13  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14  FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License along with
17  this program; if not, write to the Free Software Foundation, Inc., 59 Temple
18  Place, Suite 330, Boston, MA 02111-1307 USA.
19 
20  Turku PET Centre hereby disclaims all copyright interest in the program.
21  Juhani Knuuti
22  Director, Professor
23  Turku PET Centre, Turku, Finland, http://www.turkupetcentre.fi/
24 
25  Assumptions:
26  1. All data is always saved in little endian byte order (i386 and VAX).
27  2. Data is automatically saved in one of the little endian formats
28  as specified in header data_type.
29  3. VAX data can be saved correctly only in 2-byte formats.
30 
31  Modification history:
32  2003-07-21 Vesa Oikonen
33  Created based on previous contents of ecat63.c.
34  2003-09-08 VO
35  Added function ecat63WriteImageMatrix() and ecat63WriteScanMatrix().
36  2003-11-30 VO
37  For now, calls temp_roundf() instead of roundf().
38  2004-01-07 VO
39  ecat63WriteImageMatrix(): corrected img min&max in header.
40  2004-09-17 VO
41  Doxygen style comments.
42  2004-12-28 VO
43  Included function ecat63_is_scaling_needed().
44  This function is applied to determine whether scal factor can be set to
45  one in case that all pixel values are close to integers and small enough.
46  2007-01-27 VO
47  Unsigned char pointer was corrected to signed in ecat63WriteMatdata().
48  2007-09-02 VO
49  Backup file extension changed from % to .bak.
50 
51 ******************************************************************************/
52 #include <stdio.h>
53 #include <stdlib.h>
54 #include <math.h>
55 #include <string.h>
56 #include <unistd.h>
57 #include <time.h>
58 /*****************************************************************************/
59 #include <swap.h>
60 #include <petc99.h>
61 #include "include/ecat63.h"
62 /*****************************************************************************/
63 
64 /*****************************************************************************/
74  char buf[MatBLKSIZE];
75  int i, little, tovax;
76 
77 
78  if(ECAT63_TEST) printf("ecat63WriteMainheader()\n");
79  little=little_endian();
80  /* Clear buf */
81  memset(buf, 0, MatBLKSIZE);
82  /* Check arguments */
83  if(fp==NULL || h->data_type<1 || h->data_type>7) return(1);
84  if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
85  tovax=1; else tovax=0;
86 
87  /* Copy short ints to buf */
88  memcpy(buf+50, &h->data_type, 2); memcpy(buf+48, &h->sw_version, 2);
89  memcpy(buf+52, &h->system_type, 2); memcpy(buf+54, &h->file_type, 2);
90  memcpy(buf+66, &h->scan_start_day, 2); memcpy(buf+68, &h->scan_start_month, 2);
91  memcpy(buf+70, &h->scan_start_year, 2); memcpy(buf+72, &h->scan_start_hour, 2);
92  memcpy(buf+74, &h->scan_start_minute, 2); memcpy(buf+76, &h->scan_start_second, 2);
93  memcpy(buf+134, &h->rot_source_speed, 2); memcpy(buf+136, &h->wobble_speed, 2);
94  memcpy(buf+138, &h->transm_source_type, 2); memcpy(buf+148, &h->transaxial_samp_mode, 2);
95  memcpy(buf+150, &h->coin_samp_mode, 2); memcpy(buf+152, &h->axial_samp_mode, 2);
96  memcpy(buf+158, &h->calibration_units, 2); memcpy(buf+160, &h->compression_code, 2);
97  memcpy(buf+350, &h->acquisition_type, 2); memcpy(buf+352, &h->bed_type, 2);
98  memcpy(buf+354, &h->septa_type, 2); memcpy(buf+376, &h->num_planes, 2);
99  memcpy(buf+378, &h->num_frames, 2); memcpy(buf+380, &h->num_gates, 2);
100  memcpy(buf+382, &h->num_bed_pos, 2); memcpy(buf+452, &h->lwr_sctr_thres, 2);
101  memcpy(buf+454, &h->lwr_true_thres, 2); memcpy(buf+456, &h->upr_true_thres, 2);
102  memcpy(buf+472, h->fill2, 40);
103  /* big to little endian if necessary */
104  if(!little) swabip(buf, MatBLKSIZE);
105 
106  /* Copy floats to buf */
107  ecat63wFloat(&h->isotope_halflife, buf+86, tovax, little);
108  ecat63wFloat(&h->gantry_tilt, buf+122, tovax, little);
109  ecat63wFloat(&h->gantry_rotation, buf+126, tovax, little);
110  ecat63wFloat(&h->bed_elevation, buf+130, tovax, little);
111  ecat63wFloat(&h->axial_fov, buf+140, tovax, little);
112  ecat63wFloat(&h->transaxial_fov, buf+144, tovax, little);
113  ecat63wFloat(&h->calibration_factor, buf+154, tovax, little);
114  ecat63wFloat(&h->init_bed_position, buf+384, tovax, little);
115  for(i=0; i<15; i++) ecat63wFloat(&h->bed_offset[i], buf+388+4*i, tovax, little);
116  ecat63wFloat(&h->plane_separation, buf+448, tovax, little);
117  ecat63wFloat(&h->init_bed_position, buf+458, tovax, little);
118 
119  /* Copy chars */
120  /*memcpy(buf+0, h->ecat_format, 14);*/
121  memcpy(buf+14, h->fill1, 14);
122  memcpy(buf+28, h->original_file_name, 20); memcpy(buf+56, h->node_id, 10);
123  memcpy(buf+78, h->isotope_code, 8); memcpy(buf+90, h->radiopharmaceutical, 32);
124  memcpy(buf+162, h->study_name, 12); memcpy(buf+174, h->patient_id, 16);
125  memcpy(buf+190, h->patient_name, 32); buf[222]=h->patient_sex;
126  memcpy(buf+223, h->patient_age, 10); memcpy(buf+233, h->patient_height, 10);
127  memcpy(buf+243, h->patient_weight, 10); buf[253]=h->patient_dexterity;
128  memcpy(buf+254, h->physician_name, 32); memcpy(buf+286, h->operator_name, 32);
129  memcpy(buf+318, h->study_description, 32); memcpy(buf+356, h->facility_name, 20);
130  memcpy(buf+462, h->user_process_code, 10);
131 
132  /* Write main header */
133  fseek(fp, 0*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=0*MatBLKSIZE) return(2);
134  if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
135 
136  return(0);
137 }
138 /*****************************************************************************/
139 
140 /*****************************************************************************/
150 int ecat63WriteImageheader(FILE *fp, int block, ECAT63_imageheader *h) {
151  char buf[MatBLKSIZE];
152  int i, little, tovax;
153 
154 
155  if(ECAT63_TEST) printf("ecat63WriteImageheader(fp, %d, ih)\n", block);
156  little=little_endian();
157  /* Clear buf */
158  memset(buf, 0, MatBLKSIZE);
159  /* Check arguments */
160  if(fp==NULL || block<3 || h->data_type<1 || h->data_type>7) return(1);
161  if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
162  tovax=1; else tovax=0;
163 
164  /* Copy short ints to buf */
165  memcpy(buf+126, &h->data_type, 2); memcpy(buf+128, &h->num_dimensions, 2);
166  memcpy(buf+132, &h->dimension_1, 2); memcpy(buf+134, &h->dimension_2, 2);
167  memcpy(buf+176, &h->image_min, 2); memcpy(buf+178, &h->image_max, 2);
168  memcpy(buf+200, &h->slice_location, 2); memcpy(buf+202, &h->recon_start_hour, 2);
169  memcpy(buf+204, &h->recon_start_min, 2); memcpy(buf+206, &h->recon_start_sec, 2);
170  memcpy(buf+236, &h->filter_code, 2); memcpy(buf+376, &h->processing_code, 2);
171  memcpy(buf+380, &h->quant_units, 2); memcpy(buf+382, &h->recon_start_day, 2);
172  memcpy(buf+384, &h->recon_start_month, 2); memcpy(buf+386, &h->recon_start_year, 2);
173  memcpy(buf+460, h->fill2, 52);
174  /* big to little endian if necessary */
175  if(!little) swabip(buf, MatBLKSIZE);
176 
177  /* Copy floats to buf */
178  ecat63wFloat(&h->x_origin, buf+160, tovax, little);
179  ecat63wFloat(&h->y_origin, buf+164, tovax, little);
180  ecat63wFloat(&h->recon_scale, buf+168, tovax, little);
181  ecat63wFloat(&h->quant_scale, buf+172, tovax, little);
182  ecat63wFloat(&h->pixel_size, buf+184, tovax, little);
183  ecat63wFloat(&h->slice_width, buf+188, tovax, little);
184  ecat63wFloat(&h->image_rotation, buf+296, tovax, little);
185  ecat63wFloat(&h->plane_eff_corr_fctr, buf+300, tovax, little);
186  ecat63wFloat(&h->decay_corr_fctr, buf+304, tovax, little);
187  ecat63wFloat(&h->loss_corr_fctr, buf+308, tovax, little);
188  ecat63wFloat(&h->intrinsic_tilt, buf+312, tovax, little);
189  ecat63wFloat(&h->ecat_calibration_fctr, buf+388, tovax, little);
190  ecat63wFloat(&h->well_counter_cal_fctr, buf+392, tovax, little);
191  for(i=0; i<6; i++) ecat63wFloat(&h->filter_params[i], buf+396+4*i, tovax, little);
192 
193  /* Copy ints to buf */
194  ecat63wInt(&h->frame_duration, buf+192, tovax, little);
195  ecat63wInt(&h->frame_start_time, buf+196, tovax, little);
196  ecat63wInt(&h->scan_matrix_num, buf+238, tovax, little);
197  ecat63wInt(&h->norm_matrix_num, buf+242, tovax, little);
198  ecat63wInt(&h->atten_cor_mat_num, buf+246, tovax, little);
199 
200  /* Copy chars */
201  memcpy(buf+0, h->fill1, 126); memcpy(buf+420, h->annotation, 40);
202 
203  /* Write subheader */
204  fseek(fp, (block-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(block-1)*MatBLKSIZE) return(2);
205  if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
206 
207  return(0);
208 }
209 /*****************************************************************************/
210 
211 /*****************************************************************************/
221 int ecat63WriteAttnheader(FILE *fp, int block, ECAT63_attnheader *h) {
222  unsigned char buf[MatBLKSIZE];
223  int little, tovax;
224 
225  if(ECAT63_TEST) printf("ecat63WriteAttnheader(fp, %d, ah)\n", block);
226  little=little_endian();
227  /* Clear buf */
228  memset(buf, 0, MatBLKSIZE);
229  /* Check arguments */
230  if(fp==NULL || block<3 || h->data_type<1 || h->data_type>7) return(1);
231  if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
232  tovax=1; else tovax=0;
233 
234  /* Copy short ints to buf */
235  memcpy(buf+126, &h->data_type, 2); memcpy(buf+128, &h->attenuation_type, 2);
236  memcpy(buf+132, &h->dimension_1, 2); memcpy(buf+134, &h->dimension_2, 2);
237  /* big to little endian if necessary */
238  if(!little) swabip(buf, MatBLKSIZE);
239 
240  /* Copy floats to buf */
241  ecat63wFloat(&h->scale_factor, buf+182, tovax, little);
242  ecat63wFloat(&h->x_origin, buf+186, tovax, little);
243  ecat63wFloat(&h->y_origin, buf+190, tovax, little);
244  ecat63wFloat(&h->x_radius, buf+194, tovax, little);
245  ecat63wFloat(&h->y_radius, buf+198, tovax, little);
246  ecat63wFloat(&h->tilt_angle, buf+202, tovax, little);
247  ecat63wFloat(&h->attenuation_coeff, buf+206, tovax, little);
248  ecat63wFloat(&h->sample_distance, buf+210, tovax, little);
249 
250  /* Write subheader */
251  fseek(fp, (block-1)*MatBLKSIZE, SEEK_SET);
252  if(ftell(fp)!=(block-1)*MatBLKSIZE) return(2);
253  if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
254 
255  return(0);
256 }
257 /*****************************************************************************/
258 
259 /*****************************************************************************/
269 int ecat63WriteScanheader(FILE *fp, int block, ECAT63_scanheader *h) {
270  unsigned char buf[MatBLKSIZE];
271  int i, little, tovax;
272 
273 
274  if(ECAT63_TEST) printf("ecat63WriteScanheader(fp, %d, ih)\n", block);
275  little=little_endian();
276  /* Clear buf */
277  memset(buf, 0, MatBLKSIZE);
278  /* Check arguments */
279  if(fp==NULL || block<3 || h->data_type<1 || h->data_type>7) return(1);
280  if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
281  tovax=1; else tovax=0;
282 
283  /* Copy short ints to buf */
284  memcpy(buf+126, &h->data_type, 2);
285  memcpy(buf+132, &h->dimension_1, 2); memcpy(buf+134, &h->dimension_2, 2);
286  memcpy(buf+136, &h->smoothing, 2); memcpy(buf+138, &h->processing_code, 2);
287  memcpy(buf+170, &h->frame_duration_sec, 2);
288  memcpy(buf+192, &h->scan_min, 2); memcpy(buf+194, &h->scan_max, 2);
289  memcpy(buf+468, h->fill2, 44);
290  /* big to little endian if necessary */
291  if(!little) swabip(buf, MatBLKSIZE);
292 
293  /* Copy floats to buf */
294  ecat63wFloat(&h->sample_distance, buf+146, tovax, little);
295  ecat63wFloat(&h->isotope_halflife, buf+166, tovax, little);
296  ecat63wFloat(&h->scale_factor, buf+182, tovax, little);
297  for(i=0; i<16; i++) {
298  ecat63wFloat(&h->cor_singles[i], buf+316+4*i, tovax, little);
299  ecat63wFloat(&h->uncor_singles[i], buf+380+4*i, tovax, little);
300  }
301  ecat63wFloat(&h->tot_avg_cor, buf+444, tovax, little);
302  ecat63wFloat(&h->tot_avg_uncor, buf+448, tovax, little);
303  ecat63wFloat(&h->loss_correction_fctr, buf+464, tovax, little);
304 
305  /* Copy ints to buf */
306  ecat63wInt(&h->gate_duration, buf+172, tovax, little);
307  ecat63wInt(&h->r_wave_offset, buf+176, tovax, little);
308  ecat63wInt(&h->prompts, buf+196, tovax, little);
309  ecat63wInt(&h->delayed, buf+200, tovax, little);
310  ecat63wInt(&h->multiples, buf+204, tovax, little);
311  ecat63wInt(&h->net_trues, buf+208, tovax, little);
312  ecat63wInt(&h->total_coin_rate, buf+452, tovax, little);
313  ecat63wInt(&h->frame_start_time, buf+456, tovax, little);
314  ecat63wInt(&h->frame_duration, buf+460, tovax, little);
315 
316  /* Copy chars */
317  memcpy(buf+0, h->fill1, 126);
318 
319  /* Write subheader */
320  fseek(fp, (block-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(block-1)*MatBLKSIZE) return(2);
321  if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
322 
323  return(0);
324 }/*****************************************************************************/
325 
326 /*****************************************************************************/
336 int ecat63WriteNormheader(FILE *fp, int block, ECAT63_normheader *h)
337 {
338  unsigned char buf[MatBLKSIZE];
339  int little, tovax;
340 
341  if(ECAT63_TEST) printf("ecat63WriteNormheader(fp, %d, nh)\n", block);
342  little=little_endian();
343  /* Clear buf */
344  memset(buf, 0, MatBLKSIZE);
345  /* Check arguments */
346  if(fp==NULL || block<3 || h->data_type<1 || h->data_type>7) return(1);
347  if(h->data_type==VAX_I2 || h->data_type==VAX_I4 || h->data_type==VAX_R4)
348  tovax=1; else tovax=0;
349 
350  /* Copy short ints to buf */
351  memcpy(buf+126, &h->data_type, 2);
352  memcpy(buf+132, &h->dimension_1, 2); memcpy(buf+134, &h->dimension_2, 2);
353  memcpy(buf+372, &h->norm_hour, 2);
354  memcpy(buf+376, &h->norm_minute, 2);
355  memcpy(buf+380, &h->norm_second, 2);
356  memcpy(buf+384, &h->norm_day, 2);
357  memcpy(buf+388, &h->norm_month, 2);
358  memcpy(buf+392, &h->norm_year, 2);
359  /* big to little endian if necessary */
360  if(!little) swabip(buf, MatBLKSIZE);
361 
362  /* Copy floats to buf */
363  ecat63wFloat(&h->scale_factor, buf+182, tovax, little);
364  ecat63wFloat(&h->fov_source_width, buf+198, tovax, little);
365 
366  /* Write subheader */
367  fseek(fp, (block-1)*MatBLKSIZE, SEEK_SET);
368  if(ftell(fp)!=(block-1)*MatBLKSIZE) return(2);
369  if(fwrite(buf, 1, 1*MatBLKSIZE, fp) != 1*MatBLKSIZE) return(3);
370 
371  return(0);
372 }
373 /*****************************************************************************/
374 
375 /*****************************************************************************/
386 FILE *ecat63Create(const char *fname, ECAT63_mainheader *h) {
387  FILE *fp;
388  char tmp[FILENAME_MAX];
389  int buf[MatBLKSIZE/4];
390 
391  if(ECAT63_TEST) printf("ecat63Create()\n");
392  /* Check the arguments */
393  if(fname==NULL || h==NULL) return(NULL);
394  /* Check if file exists; backup, if necessary */
395  if(access(fname, 0) != -1) {
396  strcpy(tmp, fname); strcat(tmp, BACKUP_EXTENSION);
397  if(access(tmp, 0) != -1) remove(tmp);
398  if(ECAT63_TEST) printf("Renaming %s -> %s\n", fname, tmp);
399  rename(fname, tmp);
400  }
401  /* Open file */
402  fp=fopen(fname, "wb+"); if(fp==NULL) return(fp);
403  /* Write main header */
404  if(ecat63WriteMainheader(fp, h)) return(NULL);
405  /* Construct an empty matrix list ; convert to little endian if necessary */
406  memset(buf, 0, MatBLKSIZE);
407  buf[0]=31; buf[1]=2; if(!little_endian()) swawbip(buf, MatBLKSIZE);
408  /* Write data buffer */
409  fseek(fp, (MatFirstDirBlk-1)*MatBLKSIZE, SEEK_SET);
410  if(ftell(fp)!=(MatFirstDirBlk-1)*MatBLKSIZE) return(NULL);
411  if(fwrite(buf, 4, MatBLKSIZE/4, fp) != MatBLKSIZE/4) return(NULL);
412  /* OK, then return file pointer */
413  return(fp);
414 }
415 /*****************************************************************************/
416 
417 /*****************************************************************************/
429 int ecat63WriteImage(FILE *fp, int matnum, ECAT63_imageheader *h, void *data) {
430  int nxtblk, blkNr, data_size, pxlNr, pxlSize, ret;
431 
432  if(ECAT63_TEST) printf("ecat63WriteImage(fp, %d, ih, data)\n", matnum);
433  if(fp==NULL || matnum<1 || h==NULL || data==NULL) return(1);
434  /* nr of pixels */
435  pxlNr=h->dimension_1*h->dimension_2; if(pxlNr<1) return(2);
436  /* mem taken by one pixel */
437  switch(h->data_type) {
438  case BYTE_TYPE: pxlSize=1;
439  break;
440  case VAX_I2:
441  case SUN_I2: pxlSize=2;
442  break;
443  case VAX_I4: return(3);
444  case VAX_R4: return(3);
445  case IEEE_R4:
446  case SUN_I4: pxlSize=4;
447  break;
448  default: return(2);
449  }
450  /* mem taken by all pixels */
451  data_size=pxlNr*pxlSize;
452  /* block nr taken by all pixels */
453  blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(3);
454  /* Get block number for matrix header and data */
455  nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) return(4);
456  if(ECAT63_TEST) printf(" block=%d\n", nxtblk);
457  /* Write header */
458  ret=ecat63WriteImageheader(fp, nxtblk, h); if(ret) return(40+ret);
459  /* Write matrix data */
460  ret=ecat63WriteMatdata(fp, nxtblk+1, data, pxlNr, pxlSize);
461  if(ret) return(50+ret);
462  return 0;
463 }
464 /*****************************************************************************/
465 
466 /*****************************************************************************/
478 int ecat63WriteScan(FILE *fp, int matnum, ECAT63_scanheader *h, void *data) {
479  int nxtblk, blkNr, data_size, pxlNr, pxlSize, ret;
480 
481  if(ECAT63_TEST) printf("ecat63WriteScan(fp, %d, sh, data)\n", matnum);
482  if(fp==NULL || matnum<1 || h==NULL || data==NULL) return(1);
483  /* nr of pixels */
484  pxlNr=h->dimension_1*h->dimension_2; if(pxlNr<1) return(1);
485  /* mem taken by one pixel */
486  switch(h->data_type) {
487  case BYTE_TYPE: pxlSize=1;
488  break;
489  case VAX_I2:
490  case SUN_I2: pxlSize=2;
491  break;
492  case VAX_I4: return(3);
493  case VAX_R4: return(3);
494  case IEEE_R4:
495  case SUN_I4: pxlSize=4;
496  break;
497  default: return(2);
498  }
499  /* mem taken by all pixels */
500  data_size=pxlNr*pxlSize;
501  /* block nr taken by all pixels */
502  blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(3);
503  /* Get block number for matrix header and data */
504  nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) return(4);
505  if(ECAT63_TEST) printf(" block=%d\n", nxtblk);
506  /* Write header */
507  ret=ecat63WriteScanheader(fp, nxtblk, h); if(ret) return(40+ret);
508  /* Write matrix data */
509  ret=ecat63WriteMatdata(fp, nxtblk+1, data, pxlNr, pxlSize);
510  if(ret) return(50+ret);
511  return 0;
512 }
513 /*****************************************************************************/
514 
515 /*****************************************************************************/
527 int ecat63WriteNorm(FILE *fp, int matnum, ECAT63_normheader *h, void *data) {
528  int nxtblk, blkNr, data_size, pxlNr, pxlSize, ret;
529 
530  if(ECAT63_TEST) printf("ecat63WriteNorm(fp, %d, nh, data)\n", matnum);
531  if(fp==NULL || matnum<1 || h==NULL || data==NULL) return(1);
532  /* nr of pixels */
533  pxlNr=h->dimension_1*h->dimension_2; if(pxlNr<1) return(1);
534  /* mem taken by one pixel */
535  switch(h->data_type) {
536  case BYTE_TYPE: pxlSize=1;
537  break;
538  case VAX_I2:
539  case SUN_I2: pxlSize=2;
540  break;
541  case VAX_I4: return(3);
542  case VAX_R4: return(3);
543  case IEEE_R4:
544  case SUN_I4: pxlSize=4;
545  break;
546  default: return(2);
547  }
548  /* mem taken by all pixels */
549  data_size=pxlNr*pxlSize;
550  /* block nr taken by all pixels */
551  blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(3);
552  /* Get block number for matrix header and data */
553  nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) return(4);
554  if(ECAT63_TEST) printf(" block=%d\n", nxtblk);
555  /* Write header */
556  ret=ecat63WriteNormheader(fp, nxtblk, h); if(ret) return(40+ret);
557  /* Write matrix data */
558  ret=ecat63WriteMatdata(fp, nxtblk+1, data, pxlNr, pxlSize);
559  if(ret) return(50+ret);
560  return 0;
561 }
562 /*****************************************************************************/
563 
564 /*****************************************************************************/
576 int ecat63WriteAttn(FILE *fp, int matnum, ECAT63_attnheader *h, void *data) {
577  int nxtblk, blkNr, data_size, pxlNr, pxlSize, ret;
578 
579  if(ECAT63_TEST) printf("ecat63WriteAttn(fp, %d, ah, data)\n", matnum);
580  if(fp==NULL || matnum<1 || h==NULL || data==NULL) return(1);
581  /* nr of pixels */
582  pxlNr=h->dimension_1*h->dimension_2; if(pxlNr<1) return(1);
583  /* mem taken by one pixel */
584  switch(h->data_type) {
585  case BYTE_TYPE: pxlSize=1;
586  break;
587  case VAX_I2:
588  case SUN_I2: pxlSize=2;
589  break;
590  case VAX_I4: return(3);
591  case VAX_R4: return(3);
592  case IEEE_R4:
593  case SUN_I4: pxlSize=4;
594  break;
595  default: return(2);
596  }
597  /* mem taken by all pixels */
598  data_size=pxlNr*pxlSize;
599  /* block nr taken by all pixels */
600  blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(3);
601  /* Get block number for matrix header and data */
602  nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) return(4);
603  if(ECAT63_TEST) printf(" block=%d\n", nxtblk);
604  /* Write header */
605  ret=ecat63WriteAttnheader(fp, nxtblk, h); if(ret) return(40+ret);
606  /* Write matrix data */
607  ret=ecat63WriteMatdata(fp, nxtblk+1, data, pxlNr, pxlSize);
608  if(ret) return(50+ret);
609  return 0;
610 }
611 /*****************************************************************************/
612 
613 /*****************************************************************************/
629 int ecat63WriteMatdata(FILE *fp, int strtblk, char *data, int pxlNr, int pxlSize) {
630  unsigned char buf[MatBLKSIZE];
631  char *dptr;
632  int i, blkNr, dataSize, byteNr, little;
633 
634  if(ECAT63_TEST) printf("ecat63WriteMatdata(fp, %d, data, %d, %d)\n", strtblk, pxlNr, pxlSize);
635  if(fp==NULL || strtblk<1 || data==NULL || pxlNr<1 || pxlSize<1) return(1);
636  little=little_endian(); memset(buf, 0, MatBLKSIZE);
637  dataSize=pxlNr*pxlSize; if(dataSize<1) return(1);
638  /* block nr taken by all pixels */
639  blkNr=(dataSize+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) return(1);
640  if(ECAT63_TEST>1) printf(" blkNr=%d\n", blkNr);
641  /* Search the place for writing */
642  fseek(fp, (strtblk-1)*MatBLKSIZE, SEEK_SET); if(ftell(fp)!=(strtblk-1)*MatBLKSIZE) return(2);
643  /* Save blocks one at a time */
644  for(i=0, dptr=data; i<blkNr && dataSize>0; i++) {
645  byteNr=(dataSize<MatBLKSIZE)?dataSize:MatBLKSIZE;
646  memcpy(buf, dptr, byteNr);
647  /* Change matrix byte order in big endian platforms */
648  if(!little_endian()) {
649  if(pxlSize==2) swabip(buf, byteNr);
650  else if(pxlSize==4) swawbip(buf, byteNr);
651  }
652  /* Write block */
653  if(fwrite(buf, 1, MatBLKSIZE, fp)!=MatBLKSIZE) return(3);
654  /* Prepare for the next block */
655  dptr+=byteNr;
656  dataSize-=byteNr;
657  } /* next block */
658  return(0);
659 }
660 /*****************************************************************************/
661 
662 /*****************************************************************************/
672 int ecat63_is_scaling_needed(float amax, float *data, int nr) {
673  int i;
674  double d;
675 
676  if(nr<1 || data==NULL) return(0);
677  /* scaling is necessary if all values are between -1 - 1 */
678  if(amax<0.9999) return(1);
679  /* Lets check first if at least the max value is close to integers or not */
680  if(modf(amax, &d)>0.0001) return(1);
681  /* if it is, then check all pixels */
682  for(i=0; i<nr; i++) if(modf(*data++, &d)>0.0001) return(1);
683  return(0);
684 }
685 /*****************************************************************************/
686 
687 /*****************************************************************************/
700 int ecat63WriteImageMatrix(FILE *fp, int matnum, ECAT63_imageheader *h, float *fdata) {
701  int i, nxtblk, blkNr, data_size, pxlNr, ret;
702  float *fptr, fmin, fmax, g, f;
703  char *mdata, *mptr;
704  short int *sptr;
705 
706 
707 
708  if(ECAT63_TEST) printf("ecat63WriteImageMatrix(fp, %d, h, data)\n", matnum);
709  if(fp==NULL || matnum<1 || h==NULL || fdata==NULL) {
710  sprintf(ecat63errmsg, "invalid function parameter.\n");
711  return(1);
712  }
713  /* nr of pixels */
714  pxlNr=h->dimension_1*h->dimension_2;
715  if(pxlNr<1) {
716  sprintf(ecat63errmsg, "invalid matrix dimension.\n");
717  return(3);
718  }
719  /* How much memory is needed for ALL pixels */
720  data_size=pxlNr*ecat63pxlbytes(h->data_type);
721  /* block nr taken by all pixels */
722  blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) {
723  sprintf(ecat63errmsg, "invalid block number.\n");
724  return(4);
725  }
726  /* Allocate memory for matrix data */
727  mdata=(char*)calloc(blkNr, MatBLKSIZE); if(mdata==NULL) {
728  sprintf(ecat63errmsg, "out of memory.\n");
729  return(5);
730  }
731  /* Search for min and max for calculation of scale factor */
732  fptr=fdata; fmin=fmax=*fptr;
733  for(i=0; i<pxlNr; i++, fptr++) {
734  if(*fptr>fmax) fmax=*fptr; else if(*fptr<fmin) fmin=*fptr;
735  }
736  if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
737  if(g>0) f=32766./g; else f=1.0;
738  /* Check if pixels values can be left as such with scale_factor = 1 */
739  fptr=fdata;
740  if(f>=1.0 && ecat63_is_scaling_needed(g, fptr, pxlNr)==0) f=1.0;
741  /* Scale matrix data to shorts */
742  h->quant_scale=1.0/f;
743  sptr=(short int*)mdata; fptr=fdata;
744  for(i=0; i<pxlNr; i++, sptr++, fptr++) *sptr=(short int)temp_roundf(f*(*fptr));
745  /* Set header short min & max */
746  h->image_min=(short int)temp_roundf(f*fmin);
747  h->image_max=(short int)temp_roundf(f*fmax);
748  /* Get block number for matrix header and data */
749  nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) {
750  sprintf(ecat63errmsg, "cannot determine matrix block (%d).\n", -nxtblk);
751  free(mdata); return(8);
752  }
753  if(ECAT63_TEST>2) printf(" block=%d fmin=%g fmax=%g\n", nxtblk, fmin, fmax);
754  /* Write header */
755  ret=ecat63WriteImageheader(fp, nxtblk, h); if(ret) {
756  sprintf(ecat63errmsg, "cannot write subheader (%d).\n", ret);
757  free(mdata); return(10);
758  }
759  /* Write matrix data */
760  mptr=mdata;
761  ret=ecat63WriteMatdata(fp, nxtblk+1, mptr, pxlNr, ecat63pxlbytes(h->data_type));
762  free(mdata);
763  if(ret) {
764  sprintf(ecat63errmsg, "cannot write matrix data (%d).\n", ret);
765  return(13);
766  }
767  return(0);
768 }
769 /*****************************************************************************/
770 
771 /*****************************************************************************/
784 int ecat63WriteScanMatrix(FILE *fp, int matnum, ECAT63_scanheader *h, float *fdata) {
785  int i, nxtblk, blkNr, data_size, pxlNr, ret;
786  float *fptr, fmin, fmax, g, f;
787  char *mdata, *mptr;
788  short int *sptr;
789 
790 
791  if(ECAT63_TEST) printf("ecat63WriteScanMatrix(fp, %d, h, data)\n", matnum);
792  if(fp==NULL || matnum<1 || h==NULL || fdata==NULL) {
793  sprintf(ecat63errmsg, "invalid function parameter.\n");
794  return(1);
795  }
796  /* nr of pixels */
797  pxlNr=h->dimension_1*h->dimension_2;
798  if(pxlNr<1) {
799  sprintf(ecat63errmsg, "invalid matrix dimension.\n");
800  return(3);
801  }
802  /* How much memory is needed for ALL pixels */
803  data_size=pxlNr*ecat63pxlbytes(h->data_type);
804  /* block nr taken by all pixels */
805  blkNr=(data_size+MatBLKSIZE-1)/MatBLKSIZE; if(blkNr<1) {
806  sprintf(ecat63errmsg, "invalid block number.\n");
807  return(4);
808  }
809  /* Allocate memory for matrix data */
810  mdata=(char*)calloc(blkNr, MatBLKSIZE); if(mdata==NULL) {
811  sprintf(ecat63errmsg, "out of memory.\n");
812  return(5);
813  }
814  /* Search for min and max for calculation of scale factor */
815  fptr=fdata; fmin=fmax=*fptr;
816  for(i=0; i<pxlNr; i++, fptr++) {
817  if(*fptr>fmax) fmax=*fptr; else if(*fptr<fmin) fmin=*fptr;
818  }
819  if(fabs(fmin)>fabs(fmax)) g=fabs(fmin); else g=fabs(fmax);
820  if(g>0) f=32766./g; else f=1.0;
821  /* Check if pixels values can be left as such with scale_factor = 1 */
822  fptr=fdata;
823  if(f>=1.0 && ecat63_is_scaling_needed(g, fptr, pxlNr)==0) f=1.0;
824  /* Scale matrix data to shorts */
825  h->scale_factor=1.0/f;
826  sptr=(short int*)mdata; fptr=fdata;
827  for(i=0; i<pxlNr; i++, sptr++, fptr++) *sptr=(short int)temp_roundf(f*(*fptr));
828  /* Set header short min & max */
829  h->scan_min=(short int)temp_roundf(f*fmin);
830  h->scan_max=(short int)temp_roundf(f*fmax);
831  /* Get block number for matrix header and data */
832  nxtblk=ecat63Matenter(fp, matnum, blkNr); if(nxtblk<1) {
833  sprintf(ecat63errmsg, "cannot determine matrix block (%d).\n", -nxtblk);
834  free(mdata); return(8);
835  }
836  if(ECAT63_TEST>2) printf(" block=%d fmin=%g fmax=%g\n", nxtblk, fmin, fmax);
837  /* Write header */
838  ret=ecat63WriteScanheader(fp, nxtblk, h); if(ret) {
839  sprintf(ecat63errmsg, "cannot write subheader (%d).\n", ret);
840  free(mdata); return(10);
841  }
842  /* Write matrix data */
843  mptr=mdata;
844  ret=ecat63WriteMatdata(fp, nxtblk+1, mptr, pxlNr, ecat63pxlbytes(h->data_type));
845  free(mdata);
846  if(ret) {
847  sprintf(ecat63errmsg, "cannot write matrix data (%d).\n", ret);
848  return(13);
849  }
850  return(0);
851 }
852 /*****************************************************************************/
853 
854 /*****************************************************************************/
863 void ecat63wFloat(float *bufi, void *bufo, int tovax, int islittle) {
864  unsigned int ul;
865 
866  memcpy(&ul, bufi, 4); if(ul==0) {memcpy(bufo, bufi, 4); return;}
867  if(tovax) { /* If VAX format is needed */
868  ul+=(2L<<23); /* increase exp by 2 */
869  /* Swap words on i386 and bytes on SUN */
870  if(islittle) swawip(&ul, 4); else swabip(&ul, 4);
871  } else {
872  if(!islittle) swawbip(&ul, 4); /* Switch words and bytes on SUN */
873  }
874  memcpy(bufo, &ul, 4);
875 }
885 void ecat63wInt(int *bufi, void *bufo, int tovax, int islittle) {
886  int i;
887 
888  /* Swap both words and bytes on SUN */
889  memcpy(&i, bufi, 4); if(!islittle) swawbip(&i, 4);
890  memcpy(bufo, &i, 4);
891 }
892 /*****************************************************************************/
893 
894 /*****************************************************************************/