My Project
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/*****************************************************************************/
150int 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/*****************************************************************************/
221int 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/*****************************************************************************/
269int 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/*****************************************************************************/
336int 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/*****************************************************************************/
386FILE *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/*****************************************************************************/
429int 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/*****************************************************************************/
478int 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/*****************************************************************************/
527int 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/*****************************************************************************/
576int 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/*****************************************************************************/
629int 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/*****************************************************************************/
672int 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/*****************************************************************************/
700int 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/*****************************************************************************/
784int 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/*****************************************************************************/
863void 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}
885void 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/*****************************************************************************/
#define BACKUP_EXTENSION
Definition: analyze.h:19
#define MatFirstDirBlk
Definition: ecat63.h:28
#define SUN_I4
Definition: ecat63.h:36
char ecat63errmsg[128]
Definition: ecat63.h:50
int ECAT63_TEST
Definition: ecat63.h:52
#define VAX_R4
Definition: ecat63.h:33
#define IEEE_R4
Definition: ecat63.h:34
#define BYTE_TYPE
Definition: ecat63.h:30
#define VAX_I4
Definition: ecat63.h:32
#define VAX_I2
Definition: ecat63.h:31
#define MatBLKSIZE
Definition: ecat63.h:27
#define SUN_I2
Definition: ecat63.h:35
int ecat63Matenter(FILE *fp, int matnum, int blkNr)
Definition: ecat63ml.c:186
int ecat63pxlbytes(short int data_type)
Definition: ecat63r.c:711
int ecat63WriteScan(FILE *fp, int matnum, ECAT63_scanheader *h, void *data)
Definition: ecat63w.c:478
int ecat63WriteImageheader(FILE *fp, int block, ECAT63_imageheader *h)
Definition: ecat63w.c:150
int ecat63WriteScanMatrix(FILE *fp, int matnum, ECAT63_scanheader *h, float *fdata)
Definition: ecat63w.c:784
int ecat63WriteAttnheader(FILE *fp, int block, ECAT63_attnheader *h)
Definition: ecat63w.c:221
void ecat63wFloat(float *bufi, void *bufo, int tovax, int islittle)
Definition: ecat63w.c:863
int ecat63WriteScanheader(FILE *fp, int block, ECAT63_scanheader *h)
Definition: ecat63w.c:269
int ecat63_is_scaling_needed(float amax, float *data, int nr)
Definition: ecat63w.c:672
int ecat63WriteAttn(FILE *fp, int matnum, ECAT63_attnheader *h, void *data)
Definition: ecat63w.c:576
int ecat63WriteNorm(FILE *fp, int matnum, ECAT63_normheader *h, void *data)
Definition: ecat63w.c:527
int ecat63WriteImageMatrix(FILE *fp, int matnum, ECAT63_imageheader *h, float *fdata)
Definition: ecat63w.c:700
int ecat63WriteNormheader(FILE *fp, int block, ECAT63_normheader *h)
Definition: ecat63w.c:336
int ecat63WriteImage(FILE *fp, int matnum, ECAT63_imageheader *h, void *data)
Definition: ecat63w.c:429
FILE * ecat63Create(const char *fname, ECAT63_mainheader *h)
Definition: ecat63w.c:386
void ecat63wInt(int *bufi, void *bufo, int tovax, int islittle)
Definition: ecat63w.c:885
int ecat63WriteMatdata(FILE *fp, int strtblk, char *data, int pxlNr, int pxlSize)
Definition: ecat63w.c:629
int ecat63WriteMainheader(FILE *fp, ECAT63_mainheader *h)
Definition: ecat63w.c:73
float tilt_angle
Definition: ecat63.h:158
float attenuation_coeff
Definition: ecat63.h:159
float x_origin
Definition: ecat63.h:158
float x_radius
Definition: ecat63.h:158
float y_origin
Definition: ecat63.h:158
float scale_factor
Definition: ecat63.h:157
float sample_distance
Definition: ecat63.h:160
float y_radius
Definition: ecat63.h:158
short int data_type
Definition: ecat63.h:155
short int dimension_1
Definition: ecat63.h:156
short int dimension_2
Definition: ecat63.h:156
short int attenuation_type
Definition: ecat63.h:155
short int data_type
Definition: ecat63.h:107
short int recon_start_day
Definition: ecat63.h:120
int frame_start_time
Definition: ecat63.h:111
short int quant_units
Definition: ecat63.h:119
short int num_dimensions
Definition: ecat63.h:107
float decay_corr_fctr
Definition: ecat63.h:118
float well_counter_cal_fctr
Definition: ecat63.h:121
short int fill2[26]
Definition: ecat63.h:123
short int image_max
Definition: ecat63.h:109
float pixel_size
Definition: ecat63.h:110
float image_rotation
Definition: ecat63.h:117
char annotation[40]
Definition: ecat63.h:122
float plane_eff_corr_fctr
Definition: ecat63.h:117
short int dimension_1
Definition: ecat63.h:107
short int recon_start_month
Definition: ecat63.h:120
float recon_scale
Definition: ecat63.h:108
float ecat_calibration_fctr
Definition: ecat63.h:121
char fill1[126]
Definition: ecat63.h:106
int atten_cor_mat_num
Definition: ecat63.h:116
short int slice_location
Definition: ecat63.h:112
float loss_corr_fctr
Definition: ecat63.h:118
float slice_width
Definition: ecat63.h:110
float quant_scale
Definition: ecat63.h:108
short int dimension_2
Definition: ecat63.h:107
short int processing_code
Definition: ecat63.h:119
short int image_min
Definition: ecat63.h:109
short int filter_code
Definition: ecat63.h:115
short int recon_start_hour
Definition: ecat63.h:113
short int recon_start_sec
Definition: ecat63.h:113
float intrinsic_tilt
Definition: ecat63.h:118
short int recon_start_year
Definition: ecat63.h:120
short int recon_start_min
Definition: ecat63.h:113
float filter_params[6]
Definition: ecat63.h:121
short int num_frames
Definition: ecat63.h:97
short int upr_true_thres
Definition: ecat63.h:99
char fill1[14]
Definition: ecat63.h:73
short int lwr_sctr_thres
Definition: ecat63.h:99
short int scan_start_minute
Definition: ecat63.h:81
short int acquisition_type
Definition: ecat63.h:95
short int scan_start_hour
Definition: ecat63.h:81
short int scan_start_day
Definition: ecat63.h:80
float gantry_rotation
Definition: ecat63.h:85
char patient_age[10]
Definition: ecat63.h:92
float transaxial_fov
Definition: ecat63.h:87
char patient_weight[10]
Definition: ecat63.h:92
float gantry_tilt
Definition: ecat63.h:85
float bed_offset[15]
Definition: ecat63.h:98
short int sw_version
Definition: ecat63.h:75
char patient_sex
Definition: ecat63.h:91
short int bed_type
Definition: ecat63.h:95
short int num_gates
Definition: ecat63.h:97
char radiopharmaceutical[32]
Definition: ecat63.h:84
short int axial_samp_mode
Definition: ecat63.h:88
char node_id[10]
Definition: ecat63.h:79
char study_description[32]
Definition: ecat63.h:94
short int lwr_true_thres
Definition: ecat63.h:99
float bed_elevation
Definition: ecat63.h:85
short int data_type
Definition: ecat63.h:76
short int septa_type
Definition: ecat63.h:95
char patient_name[32]
Definition: ecat63.h:91
short int system_type
Definition: ecat63.h:77
char physician_name[32]
Definition: ecat63.h:93
float init_bed_position
Definition: ecat63.h:98
char original_file_name[20]
Definition: ecat63.h:74
short int fill2[20]
Definition: ecat63.h:102
float isotope_halflife
Definition: ecat63.h:83
char patient_dexterity
Definition: ecat63.h:93
short int transm_source_type
Definition: ecat63.h:86
short int compression_code
Definition: ecat63.h:90
float calibration_factor
Definition: ecat63.h:89
char operator_name[32]
Definition: ecat63.h:93
char patient_id[16]
Definition: ecat63.h:91
short int rot_source_speed
Definition: ecat63.h:86
short int scan_start_year
Definition: ecat63.h:80
short int wobble_speed
Definition: ecat63.h:86
short int scan_start_month
Definition: ecat63.h:80
float axial_fov
Definition: ecat63.h:87
char patient_height[10]
Definition: ecat63.h:92
char isotope_code[8]
Definition: ecat63.h:82
short int calibration_units
Definition: ecat63.h:90
short int num_bed_pos
Definition: ecat63.h:97
char study_name[12]
Definition: ecat63.h:91
short int scan_start_second
Definition: ecat63.h:81
char facility_name[20]
Definition: ecat63.h:96
char user_process_code[10]
Definition: ecat63.h:101
short int file_type
Definition: ecat63.h:78
short int num_planes
Definition: ecat63.h:97
short int coin_samp_mode
Definition: ecat63.h:88
float plane_separation
Definition: ecat63.h:98
short int transaxial_samp_mode
Definition: ecat63.h:88
short int norm_month
Definition: ecat63.h:150
short int norm_second
Definition: ecat63.h:150
short int norm_year
Definition: ecat63.h:150
short int dimension_1
Definition: ecat63.h:148
float fov_source_width
Definition: ecat63.h:151
short int norm_hour
Definition: ecat63.h:150
float scale_factor
Definition: ecat63.h:149
short int norm_minute
Definition: ecat63.h:150
short int dimension_2
Definition: ecat63.h:148
short int data_type
Definition: ecat63.h:147
short int norm_day
Definition: ecat63.h:150
float cor_singles[16]
Definition: ecat63.h:138
float uncor_singles[16]
Definition: ecat63.h:138
int frame_duration
Definition: ecat63.h:141
short int fill2[22]
Definition: ecat63.h:143
float loss_correction_fctr
Definition: ecat63.h:142
float isotope_halflife
Definition: ecat63.h:132
short int processing_code
Definition: ecat63.h:130
short int dimension_2
Definition: ecat63.h:129
short int scan_max
Definition: ecat63.h:136
int frame_start_time
Definition: ecat63.h:141
short int frame_duration_sec
Definition: ecat63.h:133
short int scan_min
Definition: ecat63.h:136
float sample_distance
Definition: ecat63.h:131
float scale_factor
Definition: ecat63.h:135
short int data_type
Definition: ecat63.h:128
int total_coin_rate
Definition: ecat63.h:140
short int dimension_1
Definition: ecat63.h:129
char fill1[126]
Definition: ecat63.h:127
short int smoothing
Definition: ecat63.h:130
float tot_avg_uncor
Definition: ecat63.h:139
float tot_avg_cor
Definition: ecat63.h:139