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 attenuation_coeff
Definition ecat63.h:159
float scale_factor
Definition ecat63.h:157
float sample_distance
Definition ecat63.h:160
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
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 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 ecat_calibration_fctr
Definition ecat63.h:121
char fill1[126]
Definition ecat63.h:106
short int slice_location
Definition ecat63.h:112
float loss_corr_fctr
Definition ecat63.h:118
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
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
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
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
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