# HG changeset patch # User bert # Date 1108578662 0 # Node ID 5bec199a517d6c890008af9874d478a6a1bdb905 # Parent 37ea9177abdaada9d37b94016d7c7d80f2a20982 Initial checkin, code from Anthonin Reilhac diff --git a/conversion/minctoecat/ecat_write.c b/conversion/minctoecat/ecat_write.c new file mode 100644 --- /dev/null +++ b/conversion/minctoecat/ecat_write.c @@ -0,0 +1,642 @@ +#include "ecat_write.h" +#include "machine_indep.h" + +#define MatBLKSIZE 512 +#define MatFirstDirBlk 2 +#define NameLen 32 +#define IDLen 16 +#define OK 0 +#define ERROR -1 +#define V7 70 + + +char* dstypecode[NumDataSetTypes] = + { "u","s","i","a","n","pm","v8","v","p8","p","i8","S","S8","N", "FS"}; + +static char* magicNumber = "MATRIX"; + +struct MatDir { + int matnum; + int strtblk; + int endblk; + int matstat; +}; + +struct matdir { + int nmats, nmax; + struct MatDir *entry; +}; + +typedef struct matdirblk { + int nfree, nextblk, prvblk, nused ; + struct MatDir matdir[31] ; +} MatDirBlk ; + + +FILE *mat_create(const char *fname,Main_header *mhead); +FILE *mat_open( const char *fname, char *fmode); +int mat_write_main_header(FILE *fptr, Main_header *header); +int map_main_header(char *bufr,Main_header *header); +int mat_wblk(FILE *fptr,int blkno, char *bufr,int nblks); +MatDirList *mat_read_directory(MatrixFile *mptr); +MatDirBlk *mat_rdirblk(MatrixFile *file, int blknum); +int write_host_data(MatrixFile *mptr, int matnum, MatrixData *data); +int matrix_find(MatrixFile *matfile, int matnum,struct MatDir *matdir); +int mat_enter(FILE *fptr, Main_header *mhptr, int matnum, int nblks); +int insert_mdir(struct MatDir matdir, MatDirList *dirlist); +int mat_write_image_subheader(FILE *fptr, Main_header *mhptr,int blknum, Image_subheader *header); +int map_image_header(char *buf, Image_subheader *header); +int matrix_freelist(MatDirList *matdirlist); +void swaw( short *from, short *to, int length); +int mat_close(FILE *fptr); +int mat_rblk(FILE *fptr, int blkno, char *bufr, int nblks); + +MatrixFile *matrix_create(const char *fname, Main_header * proto_mhptr) { + MatrixFile *mptr = NULL; + FILE *fptr, *mat_create(); + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + + fptr = mat_create(fname, proto_mhptr); + if (!fptr) return( NULL ); + + mptr = (MatrixFile *) calloc(1, sizeof(MatrixFile)); + if (!mptr) { + fclose( fptr ); + return( NULL ); + } + + mptr->fptr = fptr; + + mptr->fname = (char *) malloc(strlen(fname) + 1); + if (!mptr->fname) { + free( mptr ); + fclose( fptr ); + return( NULL ); + } + + strcpy(mptr->fname, fname); + mptr->mhptr = (Main_header *) malloc(sizeof(Main_header)); + if (!mptr->mhptr) { + free( mptr->fname ); + free( mptr ); + fclose( fptr ); + return( NULL ); + } + + memcpy(mptr->mhptr, proto_mhptr, sizeof(Main_header)); + mptr->dirlist = mat_read_directory(mptr); + if (!mptr->dirlist) { + free( mptr->fname ); + free( mptr->mhptr ); + free( mptr ); + fclose( fptr ); + return( NULL ); + } + return mptr; +} + +FILE *mat_create(const char *fname,Main_header *mhead) { + FILE *fptr; + int bufr[MatBLKSIZE/sizeof(int)]; + int ret; + + fptr = mat_open( fname, "w+"); + if (!fptr) return( NULL ); + ret = mat_write_main_header( fptr, mhead ); + if( ret != 0 ) { + mat_close( fptr); + return( NULL ); + } + memset(bufr,0,MatBLKSIZE); + bufr[0] = 31; + bufr[1] = 2; + + ret = write_matrix_data(fptr,MatFirstDirBlk,1,(char*)bufr,SunLong); + + if( ret != 0 ) { + mat_close( fptr); + return( NULL ); + } + return (fptr); +} + +FILE *mat_open( const char *fname, char *fmode) { + FILE *fopen(), *fptr; + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + fptr = fopen(fname, fmode); + return (fptr); +} + +int mat_write_main_header(FILE *fptr, Main_header *header) { + char bufr[MatBLKSIZE]; + + map_main_header(bufr, header); + return mat_wblk(fptr, 1, bufr, 1); /* write main header at block 1 */ +} + +int map_main_header(char *bufr,Main_header *header) { + + int i = 0, j = 0; + char mn[20]; + /* set magic number */ + sprintf(mn,"%s%d%s", magicNumber,header->sw_version, + dstypecode[header->file_type]); + bufWrite(mn, bufr, &i, 14); + + /* copy buffer into struct */ + bufWrite(header->original_file_name, bufr, &i, NameLen); + bufWrite_s(header->sw_version, bufr, &i); + bufWrite_s(header->system_type, bufr, &i); + bufWrite_s(header->file_type, bufr, &i); + bufWrite(header->serial_number, bufr, &i, 10); + bufWrite_u(header->scan_start_time, bufr, &i); + bufWrite(header->isotope_code, bufr, &i, 8); + bufWrite_f(header->isotope_halflife, bufr, &i); + bufWrite(header->radiopharmaceutical, bufr, &i, NameLen); + bufWrite_f(header->gantry_tilt, bufr, &i); + bufWrite_f(header->gantry_rotation, bufr, &i); + bufWrite_f(header->bed_elevation, bufr, &i); + bufWrite_f(header->intrinsic_tilt, bufr, &i); + bufWrite_s(header->wobble_speed, bufr, &i); + bufWrite_s(header->transm_source_type, bufr, &i); + bufWrite_f(header->distance_scanned, bufr, &i); + bufWrite_f(header->transaxial_fov, bufr, &i); + bufWrite_s(header->angular_compression, bufr, &i); + bufWrite_s(header->coin_samp_mode, bufr, &i); + bufWrite_s(header->axial_samp_mode, bufr, &i); + bufWrite_f(header->calibration_factor, bufr, &i); + bufWrite_s(header->calibration_units, bufr, &i); + bufWrite_s(header->calibration_units_label, bufr, &i); + bufWrite_s(header->compression_code, bufr, &i); + bufWrite(header->study_name, bufr, &i, 12); + bufWrite(header->patient_id, bufr, &i, IDLen); + bufWrite(header->patient_name, bufr, &i, NameLen); + bufWrite(header->patient_sex, bufr, &i, 1); + bufWrite(header->patient_dexterity, bufr, &i, 1); + bufWrite_f(header->patient_age, bufr, &i); + bufWrite_f(header->patient_height, bufr, &i); + bufWrite_f(header->patient_weight, bufr, &i); + bufWrite_i(header->patient_birth_date, bufr, &i); + bufWrite(header->physician_name, bufr, &i, NameLen); + bufWrite(header->operator_name, bufr, &i, NameLen); + bufWrite(header->study_description, bufr, &i, NameLen); + bufWrite_s(header->acquisition_type, bufr, &i); + bufWrite_s(header->patient_orientation, bufr, &i); + bufWrite(header->facility_name, bufr, &i, 20); + bufWrite_s(header->num_planes, bufr, &i); + bufWrite_s(header->num_frames, bufr, &i); + bufWrite_s(header->num_gates, bufr, &i); + bufWrite_s(header->num_bed_pos, bufr, &i); + bufWrite_f(header->init_bed_position, bufr, &i); + for(j = 0; j < 15; j++) + bufWrite_f(header->bed_offset[j], bufr, &i); + bufWrite_f(header->plane_separation, bufr, &i); + bufWrite_s(header->lwr_sctr_thres, bufr, &i); + bufWrite_s(header->lwr_true_thres, bufr, &i); + bufWrite_s(header->upr_true_thres, bufr, &i); + bufWrite(header->user_process_code, bufr, &i, 10); + bufWrite_s(header->acquisition_mode, bufr, &i); + bufWrite_f(header->bin_size, bufr, &i); + bufWrite_f(header->branching_fraction, bufr, &i); + bufWrite_u(header->dose_start_time, bufr, &i); + bufWrite_f(header->dosage, bufr, &i); + bufWrite_f(header->well_counter_factor, bufr, &i); + bufWrite(header->data_units, bufr, &i, 32); + bufWrite_s(header->septa_state, bufr, &i); + return 1; +} + +int mat_wblk(FILE *fptr,int blkno, char *bufr,int nblks) { + int err; + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + + /* seek to position in file */ + err = fseek(fptr, (blkno - 1) * MatBLKSIZE, 0); + if (err) return (ERROR); + + /* write matrix data */ + err = fwrite(bufr, 1, nblks * MatBLKSIZE, fptr); + if (err == -1) return (ERROR); + if (err != nblks * MatBLKSIZE) { + matrix_errno = MAT_WRITE_ERROR; + return (ERROR); + } + return (0); +} + +MatDirList *mat_read_directory(MatrixFile *mptr) { + struct MatDir matdir; + MatDirList *dirlist; + MatDirBlk *matdirblk; + int i, blknum; + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + + dirlist = (MatDirList *) calloc(1, sizeof(MatDirList)); + if (dirlist == NULL) return (NULL); + + blknum = MatFirstDirBlk; + do { + matdirblk = mat_rdirblk(mptr, blknum); + if (matdirblk == NULL) { + free(dirlist); + return (NULL); + } + for (i = 0; i < matdirblk->nused; i++) { + matdir.matnum = matdirblk->matdir[i].matnum; + matdir.strtblk = matdirblk->matdir[i].strtblk; + matdir.endblk = matdirblk->matdir[i].endblk; + matdir.matstat = matdirblk->matdir[i].matstat; + insert_mdir(matdir, dirlist); + } + blknum = matdirblk->nextblk; + free(matdirblk); + } + while (blknum != MatFirstDirBlk); + return (dirlist); +} + +MatDirBlk *mat_rdirblk(MatrixFile *file, int blknum) { + MatDirBlk *matdirblk; + int i, j, err, ndirs; + int dirbufr[MatBLKSIZE / 4]; + FILE *fptr = file->fptr; + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + + matdirblk = (MatDirBlk *) malloc(MatBLKSIZE); + if (matdirblk == NULL) + return (NULL); + + err = read_matrix_data(fptr, blknum, 1, (char *) dirbufr, SunLong); + + if (err == ERROR) { + free(matdirblk); + return (NULL); + } + matdirblk->nfree = dirbufr[0]; + matdirblk->nextblk = dirbufr[1]; + matdirblk->prvblk = dirbufr[2]; + matdirblk->nused = dirbufr[3]; + + if (matdirblk->nused > 31) { + matrix_errno = MAT_INVALID_DIRBLK; + free(matdirblk); + return (NULL); + } + ndirs = (MatBLKSIZE / 4 - 4) / 4; + for (i = 0; i < ndirs; i++) { + matdirblk->matdir[i].matnum = 0; + matdirblk->matdir[i].strtblk = 0; + matdirblk->matdir[i].endblk = 0; + matdirblk->matdir[i].matstat = 0; + } + + for (i = 0; i < matdirblk->nused; i++) { + j = i + 1; + matdirblk->matdir[i].matnum = dirbufr[j * 4 + 0]; + matdirblk->matdir[i].strtblk = dirbufr[j * 4 + 1]; + matdirblk->matdir[i].endblk = dirbufr[j * 4 + 2]; + matdirblk->matdir[i].matstat = dirbufr[j * 4 + 3]; + } + return (matdirblk); +} + +int matrix_write(MatrixFile *mptr, int matnum,MatrixData *data) { + int slice ; + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + if (mptr == NULL) matrix_errno = MAT_READ_FROM_NILFPTR ; + else if (mptr->mhptr == NULL) matrix_errno = MAT_NOMHD_FILE_OBJECT ; + else if (data->shptr == NULL) matrix_errno = MAT_NIL_SHPTR ; + else if (data->data_ptr == NULL) matrix_errno = MAT_NIL_DATA_PTR ; + if (matrix_errno != OK) return (ERROR) ; + return write_host_data(mptr, matnum, data); +} + +int write_host_data(MatrixFile *mptr, int matnum, MatrixData *data) { + struct MatDir matdir, dir_entry ; + Image_subheader *imagesub ; + int status, blkno, nblks ; + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + status = OK ; + nblks = (data->data_size+511)/512; + + if (matrix_find(mptr, matnum, &matdir) == ERROR) { + blkno = mat_enter(mptr->fptr, mptr->mhptr, matnum, nblks) ; + if( blkno == ERROR ) return( ERROR ); + dir_entry.matnum = matnum ; + dir_entry.strtblk = blkno ; + + /*le 24 Avril 2001, j'ai enleve le -1*/ + dir_entry.endblk = dir_entry.strtblk + nblks /*- 1*/ ; + dir_entry.matstat = 1 ; + insert_mdir(dir_entry, mptr->dirlist) ; + matdir = dir_entry ; + } + + imagesub = (Image_subheader *) data->shptr ; + if (imagesub == NULL) { + imagesub = (Image_subheader *) calloc(1, MatBLKSIZE); + data->shptr = (caddr_t)imagesub; + } /* use MatrixData info */ + imagesub->x_pixel_size = data->pixel_size; + imagesub->y_pixel_size = data->y_size; + imagesub->z_pixel_size = data->z_size; + imagesub->num_dimensions = 3; + imagesub->x_dimension = data->xdim; + imagesub->y_dimension = data->ydim; + imagesub->z_dimension = data->zdim; + imagesub->image_max = (int)(data->data_max/data->scale_factor); + imagesub->image_min = (int)(data->data_min/data->scale_factor); + imagesub->scale_factor = data->scale_factor; + imagesub->data_type = data->data_type; + if( mat_write_image_subheader(mptr->fptr,mptr->mhptr,matdir.strtblk, imagesub) == ERROR ) return( ERROR ); + status = write_matrix_data(mptr->fptr, matdir.strtblk+1, nblks, data->data_ptr, imagesub->data_type) ; + if( status == ERROR ) return( ERROR ); + return(status) ; +} + +int matrix_find(MatrixFile *matfile, int matnum,struct MatDir *matdir) { + MatDirNode *node ; + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + if (matfile == NULL) return(ERROR) ; + if (matfile->dirlist == NULL) return(ERROR) ; + node = matfile->dirlist->first ; + while (node != NULL) + { + if (node->matnum == matnum) + { + matdir->matnum = node->matnum ; + matdir->strtblk = node->strtblk ; + matdir->endblk = node->endblk ; + matdir->matstat = node->matstat ; + break ; + } + node = node->next ; + } + if (node != NULL) return(OK) ; + else return(ERROR) ; +} + +int mat_enter(FILE *fptr, Main_header *mhptr, int matnum, int nblks) { + int dirblk, dirbufr[128], i, nxtblk, busy, oldsize; + short sw_version = mhptr->sw_version; + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + dirblk = MatFirstDirBlk; + if( fseek(fptr, 0, 0) ) return( ERROR ); + /* + * nfs locks are very time consuming lockf( fileno(fptr), F_LOCK, 0); + */ + if (read_matrix_data(fptr, dirblk, 1, (char *) dirbufr, SunLong) == ERROR) return (ERROR); + + busy = 1; + while (busy) { + nxtblk = dirblk + 1; + for (i = 4; i < 128; i += 4) { + if (dirbufr[i] == 0) { + busy = 0; + break; + } else if (dirbufr[i] == matnum) { + oldsize = dirbufr[i + 2] - dirbufr[i + 1] + 1; + if (oldsize < nblks) { + dirbufr[i] = 0xFFFFFFFF; + write_matrix_data(fptr, dirblk, 1, (char*)dirbufr, SunLong); + nxtblk = dirbufr[i + 2] + 1; + } else { + nxtblk = dirbufr[i + 1]; + dirbufr[0]++; + dirbufr[3]--; + busy = 0; + break; + } + } else { + nxtblk = dirbufr[i + 2] + 1; + } + } + if (!busy) break; + if (dirbufr[1] != MatFirstDirBlk) { + dirblk = dirbufr[1]; + if (read_matrix_data(fptr, dirblk, 1, (char *) dirbufr, SunLong) == ERROR) return (ERROR); + } else { + dirbufr[1] = nxtblk; + if (write_matrix_data(fptr, dirblk, 1, (char *) dirbufr, SunLong) == ERROR) return (ERROR); + + dirbufr[0] = 31; + dirbufr[1] = MatFirstDirBlk; + dirbufr[2] = dirblk; + dirbufr[3] = 0; + dirblk = nxtblk; + for (i = 4; i < 128; i++) + dirbufr[i] = 0; + } + } + dirbufr[i] = matnum; + dirbufr[i + 1] = nxtblk; + dirbufr[i + 2] = nxtblk + nblks; + dirbufr[i + 3] = 1; + dirbufr[0]--; + dirbufr[3]++; + if (write_matrix_data(fptr, dirblk, 1, (char*)dirbufr, SunLong) == ERROR) return (ERROR); + + if( fseek(fptr, 0, 0) ) return( ERROR ); + /* + * nfs locks are very time consuming lockf( fileno(fptr), F_UNLOCK, 0); + */ + return (nxtblk); +} + +int insert_mdir(struct MatDir matdir, MatDirList *dirlist) { + MatDirNode *node ; + + matrix_errno = MAT_OK; + matrix_errtxt[0] = '\0'; + if (dirlist == NULL) + { + dirlist = (MatDirList *) malloc(sizeof(MatDirList)) ; + if (dirlist == NULL) return(ERROR) ; + dirlist->nmats = 0 ; + dirlist->first = NULL ; + dirlist->last = NULL ; + } + node = (MatDirNode *) malloc(sizeof(MatDirNode)) ; + if (node == NULL) return(ERROR) ; + + node->matnum = matdir.matnum ; + node->strtblk = matdir.strtblk ; + node->endblk = matdir.endblk ; + node->matstat = matdir.matstat; + node->next = NULL ; + + if (dirlist->first == NULL) /* if list was empty, add first node */ + { + dirlist->first = node ; + dirlist->last = node ; + dirlist->nmats = 1 ; + } + else + { + (dirlist->last)->next = node ; + dirlist->last = node ; + ++(dirlist->nmats) ; + } + return OK; +} + +int mat_write_image_subheader(FILE *fptr, Main_header *mhptr,int blknum, Image_subheader *header) { + char buf[MatBLKSIZE]; + map_image_header(buf, header); + return mat_wblk(fptr, blknum, buf, 1); +} + +int map_image_header(char *buf, Image_subheader *header) { + int i = 0; + bufWrite_s(header->data_type, buf, &i); + bufWrite_s(header->num_dimensions, buf, &i); + bufWrite_s(header->x_dimension, buf, &i); + bufWrite_s(header->y_dimension, buf, &i); + bufWrite_s(header->z_dimension, buf, &i); + bufWrite_f(header->z_offset, buf, &i); + bufWrite_f(header->x_offset, buf, &i); + bufWrite_f(header->y_offset, buf, &i); + bufWrite_f(header->recon_zoom, buf, &i); + bufWrite_f(header->scale_factor, buf, &i); + bufWrite_s(header->image_min, buf, &i); + bufWrite_s(header->image_max, buf, &i); + bufWrite_f(header->x_pixel_size, buf, &i); + bufWrite_f(header->y_pixel_size, buf, &i); + bufWrite_f(header->z_pixel_size, buf, &i); + bufWrite_u(header->frame_duration, buf, &i); + bufWrite_u(header->frame_start_time, buf, &i); + bufWrite_s(header->filter_code, buf, &i); + bufWrite_f(header->x_resolution, buf, &i); + bufWrite_f(header->y_resolution, buf, &i); + bufWrite_f(header->z_resolution, buf, &i); + bufWrite_f(header->num_r_elements, buf, &i); + bufWrite_f(header->num_angles, buf, &i); + bufWrite_f(header->z_rotation_angle, buf, &i); + bufWrite_f(header->decay_corr_fctr, buf, &i); + bufWrite_i(header->processing_code, buf, &i); + bufWrite_u(header->gate_duration, buf, &i); + bufWrite_i(header->r_wave_offset, buf, &i); + bufWrite_i(header->num_accepted_beats, buf, &i); + bufWrite_f(header->filter_cutoff_frequency, buf, &i); + bufWrite_f(header->filter_resolution, buf, &i); + bufWrite_f(header->filter_ramp_slope, buf, &i); + bufWrite_s(header->filter_order, buf, &i); + bufWrite_f(header->filter_scatter_fraction, buf, &i); + bufWrite_f(header->filter_scatter_slope, buf, &i); + bufWrite(header->annotation, buf, &i, 40); + bufWrite_f(header->mt_1_1, buf, &i); + bufWrite_f(header->mt_1_2, buf, &i); + bufWrite_f(header->mt_1_3, buf, &i); + bufWrite_f(header->mt_2_1, buf, &i); + bufWrite_f(header->mt_2_2, buf, &i); + bufWrite_f(header->mt_2_3, buf, &i); + bufWrite_f(header->mt_3_1, buf, &i); + bufWrite_f(header->mt_3_2, buf, &i); + bufWrite_f(header->mt_3_3, buf, &i); + bufWrite_f(header->rfilter_cutoff, buf, &i); + bufWrite_f(header->rfilter_resolution, buf, &i); + bufWrite_s(header->rfilter_code, buf, &i); + bufWrite_s(header->rfilter_order, buf, &i); + bufWrite_f(header->zfilter_cutoff, buf, &i); + bufWrite_f(header->zfilter_resolution, buf, &i); + bufWrite_s(header->zfilter_code, buf, &i); + bufWrite_s(header->zfilter_order, buf, &i); + bufWrite_f(header->mt_1_4, buf, &i); + bufWrite_f(header->mt_2_4, buf, &i); + bufWrite_f(header->mt_3_4, buf, &i); + bufWrite_s(header->scatter_type, buf, &i); + bufWrite_s(header->recon_type, buf, &i); + bufWrite_s(header->recon_views, buf, &i); + return 1; +} + +int matrix_close(MatrixFile *mptr) { + int status = OK; + matrix_errno = MAT_OK; + if (mptr->fname) strcpy(matrix_errtxt,mptr->fname); + else matrix_errtxt[0] = '\0'; + if (mptr == NULL) return status; + if (mptr->mhptr != NULL) free(mptr->mhptr) ; + if (mptr->dirlist != NULL) matrix_freelist(mptr->dirlist) ; + if (mptr->fptr) status = fclose(mptr->fptr); + if (mptr->fname) free(mptr->fname); + free(mptr); + return status; +} + +int matrix_freelist(MatDirList *matdirlist) { + MatDirNode *node, *next ; + + if (matdirlist == NULL) return OK; + if (matdirlist->first != NULL) { + node = matdirlist->first ; + do { + next = node->next ; + free(node) ; + node = next ; + } + while(next != NULL) ; + } + free(matdirlist) ; + return OK; +} + + +int mat_numcod(int frame, int plane, int gate, int data, int bed) { + return ((frame)|((bed&0xF)<<12)|((plane&0xFF)<<16)|(((plane&0x300)>>8)<<9)| + ((gate&0x3F)<<24)|((data&0x3)<<30)|((data&0x4)<<9)); +} + +void swaw( short *from, short *to, int length) { + short int temp; + int i; + + for (i=0;i check if byte count less than (nblks-1) (M. Sibomana 23-oct-1997) */ + else if( err < (nblks-1)*MatBLKSIZE ) { + matrix_errno = MAT_READ_ERROR; + return( ERROR ); + } + return( 0 ); +} diff --git a/conversion/minctoecat/ecat_write.h b/conversion/minctoecat/ecat_write.h new file mode 100644 --- /dev/null +++ b/conversion/minctoecat/ecat_write.h @@ -0,0 +1,227 @@ +#include +#include + +typedef enum { ECAT6, ECAT7, Interfile } FileFormat; + +typedef enum { + MAT_OK, + MAT_READ_ERROR, + MAT_WRITE_ERROR, + MAT_INVALID_DIRBLK, + MAT_ACS_FILE_NOT_FOUND, + MAT_INTERFILE_OPEN_ERR, + MAT_FILE_TYPE_NOT_MATCH, + MAT_READ_FROM_NILFPTR, + MAT_NOMHD_FILE_OBJECT, + MAT_NIL_SHPTR, + MAT_NIL_DATA_PTR, + MAT_MATRIX_NOT_FOUND, + MAT_UNKNOWN_FILE_TYPE, + MAT_ACS_CREATE_ERR, + MAT_BAD_ATTRIBUTE, + MAT_BAD_FILE_ACCESS_MODE, + MAT_INVALID_DIMENSION, + MAT_NO_SLICES_FOUND, + MAT_INVALID_DATA_TYPE, + MAT_INVALID_MBED_POSITION +} MatrixErrorCode; + +typedef enum { + NoData, Sinogram, PetImage, AttenCor, Normalization, + PolarMap, ByteVolume, PetVolume, ByteProjection, + PetProjection, ByteImage, Short3dSinogram, Byte3dSinogram, Norm3d, + Float3dSinogram,InterfileImage, NumDataSetTypes +} DataSetType; + +typedef enum { + UnknownMatDataType, ByteData, VAX_Ix2, VAX_Ix4, + VAX_Rx4, IeeeFloat, SunShort, SunLong, NumMatrixDataTypes, ColorData, + BitData +} MatrixDataType; + +MatrixErrorCode matrix_errno; + +char matrix_errtxt[132]; + +typedef struct XMAIN_HEAD { + char magic_number[14]; + char original_file_name[32]; + short sw_version; + short system_type; + short file_type; + char serial_number[10]; + short align_0; /* 4 byte alignment purpose */ + unsigned int scan_start_time; + char isotope_code[8]; + float isotope_halflife; + char radiopharmaceutical[32]; + float gantry_tilt; + float gantry_rotation; + float bed_elevation; + float intrinsic_tilt; + short wobble_speed; + short transm_source_type; + float distance_scanned; + float transaxial_fov; + short angular_compression; + short coin_samp_mode; + short axial_samp_mode; + short align_1; + float calibration_factor; + short calibration_units; + short calibration_units_label; + short compression_code; + char study_name[12]; + char patient_id[16]; + char patient_name[32]; + char patient_sex[1]; + char patient_dexterity[1]; + float patient_age; + float patient_height; + float patient_weight; + int patient_birth_date; + char physician_name[32]; + char operator_name[32]; + char study_description[32]; + short acquisition_type; + short patient_orientation; + char facility_name[20]; + short num_planes; + short num_frames; + short num_gates; + short num_bed_pos; + float init_bed_position; + float bed_offset[15]; + float plane_separation; + short lwr_sctr_thres; + short lwr_true_thres; + short upr_true_thres; + char user_process_code[10]; + short acquisition_mode; + short align_2; + float bin_size; + float branching_fraction; + unsigned int dose_start_time; + float dosage; + float well_counter_factor; + char data_units[32]; + short septa_state; + short align_3; +} Main_header; + +typedef struct XIMAGE_SUB { + short data_type; + short num_dimensions; + short x_dimension; + short y_dimension; + short z_dimension; + short align_0; + float z_offset; + float x_offset; + float y_offset; + float recon_zoom; + float scale_factor; + short image_min; + short image_max; + float x_pixel_size; + float y_pixel_size; + float z_pixel_size; + unsigned int frame_duration; + unsigned int frame_start_time; + short filter_code; + short align_1; + float x_resolution; + float y_resolution; + float z_resolution; + float num_r_elements; + float num_angles; + float z_rotation_angle; + float decay_corr_fctr; + int processing_code; + unsigned int gate_duration; + int r_wave_offset; + int num_accepted_beats; + float filter_cutoff_frequency; + float filter_resolution; + float filter_ramp_slope; + short filter_order; + short align_2; + float filter_scatter_fraction; + float filter_scatter_slope; + char annotation[40]; + float mt_1_1; + float mt_1_2; + float mt_1_3; + float mt_2_1; + float mt_2_2; + float mt_2_3; + float mt_3_1; + float mt_3_2; + float mt_3_3; + float rfilter_cutoff; + float rfilter_resolution; + short rfilter_code; + short rfilter_order; + float zfilter_cutoff; + float zfilter_resolution; + short zfilter_code; + short zfilter_order; + float mt_1_4; + float mt_2_4; + float mt_3_4; + short scatter_type; + short recon_type; + short recon_views; + short align_3; +} Image_subheader; + +typedef struct matdirnode { + int matnum ; + int strtblk ; + int endblk ; + int matstat ; + struct matdirnode *next ; +} MatDirNode ; + +typedef struct matdirlist { + int nmats ; + MatDirNode *first ; + MatDirNode *last ; +} MatDirList ; + +typedef struct matrix_file { + char *fname ; + Main_header *mhptr ; + MatDirList *dirlist ; + FILE *fptr ; + int acs ; + FileFormat file_format; + char **interfile_header; +} MatrixFile; + +typedef struct matrixdata { + int matnum ; /* matrix number */ + MatrixFile *matfile ; /* pointer to parent */ + DataSetType mat_type ; /* type of matrix? */ + MatrixDataType data_type ; /* type of data */ + caddr_t shptr ; /* pointer to sub-header */ + caddr_t data_ptr ; /* pointer to data */ + int data_size ; /* size of data in bytes */ + int xdim; /* dimensions of data */ + int ydim; /* y dimension */ + int zdim; /* for volumes */ + float scale_factor ; /* valid if data is int? */ + float pixel_size; /* xdim data spacing (cm) */ + float y_size; /* ydim data spacing (cm) */ + float z_size; /* zdim data spacing (cm) */ + float data_min; /* min value of data */ + float data_max; /* max value of data */ + float x_origin; /* x origin of data */ + float y_origin; /* y origin of data */ + float z_origin; /* z origin of data */ +} MatrixData ; + +MatrixFile *matrix_create(const char *fname, Main_header * proto_mhptr); +int matrix_write(MatrixFile *mptr, int matnum, MatrixData *data); +int matrix_close(MatrixFile *mptr); +int mat_numcod(int frame, int plane, int gate, int data, int bed); diff --git a/conversion/minctoecat/machine_indep.c b/conversion/minctoecat/machine_indep.c new file mode 100644 --- /dev/null +++ b/conversion/minctoecat/machine_indep.c @@ -0,0 +1,415 @@ +#include +#include +#include "ecat_write.h" + +#define OK 0 +#define ERROR -1 + +#ifdef _WIN32 +u_short ntohs(u_short us) { + u_char *p = (u_char*)&us; + return ((u_short)p[1] + (p[0] >> 8)); +} +unsigned long ntohl(unsigned long ul) { + u_char *p = (u_char*)&ul; + return ((unsigned long)p[3] + (p[2] >> 8) + (p[1] >> 16) + (p[0] >> 24)); +} +#endif + + +#if defined(__STDC__) +void SWAB(const void *from, void *to, int length) +#else +void SWAB( from, to, length) +char *from, *to; +int length; +#endif +{ + if (ntohs(1) == 1) swab(from, to, length); + else memcpy(to,from,length); +} + +void SWAW ( from, to, length) +#if defined(__STDC__) +const short *from; +#else +short *from; +#endif +short *to; +int length; +{ + if (ntohs(1) == 1) swaw((short*)from,to,length); + else memcpy(to,from,length*2); +} + +float vaxftohf( bufr, off) +unsigned short *bufr; +int off; +{ + unsigned int sign_exp, high, low, mantissa, ret; + unsigned u = (bufr[off+1] << 16) + bufr[off]; + + if (u == 0) return 0.0; + sign_exp = u & 0x0000ff80; + sign_exp = (sign_exp - 0x0100) & 0xff80; + high = u & 0x0000007f; + low = u & 0xffff0000; + mantissa = (high << 16) + (low >> 16); + sign_exp = sign_exp << 16; + ret = sign_exp + mantissa; + return *(float*)(&ret); +} + +#if defined(__alpha) || defined(_WIN32) /* LITTLE_ENDIAN : alpha, intel */ +ftovaxf(f, bufr) +float f; +unsigned short *bufr; +{ + unsigned *u, sign_exp, low, high, mantissa, ret; + float f_tmp = f; + + u = (unsigned*)(&f_tmp); + if (u == 0) { + bufr[0] = bufr[1] = 0; + return; + } + sign_exp = *u & 0xff800000; + sign_exp = sign_exp >> 16; + sign_exp = (sign_exp + 0x0110) & 0xff80; + low = *u & 0x0000ffff; + high = *u & 0x007f0000; + mantissa = (high >> 16) + (low << 16); + ret = sign_exp + mantissa; + bufr[0] = ret; + bufr[1] = ret >>16; +} +#else /* BIG ENDIAN : sun hp sgi*/ +ftovaxf(orig,number) + unsigned short number[2]; + float orig; +{ + + /* convert from sun float to vax float */ + + union { + unsigned short t[2]; + float t4; + } test; + unsigned short int exp; + + number[0] = 0; + number[1] = 0; + + test.t4 = orig; + if (test.t4 == 0.0) return; + + number[1] = test.t[1]; + + exp = ((test.t[0] & 0x7f00) + 0x0100) & 0x7f00; + test.t[0] = (test.t[0] & 0x80ff) + exp; + + number[0] = test.t[0]; + +} +#endif /* LITTLE VS BIG ENDIAN*/ + +int file_data_to_host(dptr, nblks, dtype) +char *dptr; +int nblks, dtype; +{ + int i, j; + char *tmp = NULL; + + + matrix_errno = 0; + matrix_errtxt[0] = '\0'; + if ((tmp = malloc(512)) == NULL) return ERROR; + switch(dtype) + { + case ByteData: + break; + case VAX_Ix2: + if (ntohs(1) == 1) + for (i=0, j=0; i + + +#if defined(__STDC__) || defined(__cplusplus) +#if defined(__cplusplus) +extern "C" { +#endif +void SWAB(const void *from, void *to, int length); +int file_data_to_host(char *dptr, int nblks, int dtype); +int read_matrix_data( FILE *fptr, int strtblk, int nblks, + char *dptr, int dtype); +int write_matrix_data( FILE *fptr, int strtblk, int nblks, + char *dptr, int dtype); +void bufWrite(char* s, char* buf, int* i, int len); +void bufWrite_s(short val, char* buf, int* i); +void bufWrite_i(int val, char* buf, int* i); +void bufWrite_u(unsigned int val, char* buf, int* i); +void bufWrite_f(float val, char* buf, int* i); +void bufRead(char* s, char* buf, int* i, int len); +void bufRead_s(short*, char* buf, int* i); +void bufRead_i(int*, char* buf, int* i); +void bufRead_u(unsigned int*, char* buf, int* i); +void bufRead_f(float*, char* buf, int* i); +#if defined(__cplusplus) +} +#endif /* __cpluplus */ +#else /* __STDC__ */ +extern void bufWrite(), bufWrite_s(), bufWrite_i(), bufWrite_f(); +extern void bufRead(), bufRead_s(), bufRead_i(), bufRead_f(); +#endif /* __STDC__ */ + +#endif /* machine_indep_h */ diff --git a/conversion/minctoecat/minctoecat.c b/conversion/minctoecat/minctoecat.c new file mode 100644 --- /dev/null +++ b/conversion/minctoecat/minctoecat.c @@ -0,0 +1,1007 @@ +/* ----------------------------- MNI Header ----------------------------------- +@NAME : minctoecat +@INPUT : argc, argv - command line arguments +@OUTPUT : (none) +@RETURNS : error status +@DESCRIPTION: Converts a minc format file to a CTI ECAT file. +@METHOD : +@GLOBALS : +@CALLS : +@CREATED : January 31, 2005 (Anthonin Reilhac) +@MODIFIED : +@COPYRIGHT : + Copyright 2005 Anthonin Reilhac, McConnell Brain Imaging Centre, + Montreal Neurological Institute, McGill University. + Permission to use, copy, modify, and distribute this + software and its documentation for any purpose and without + fee is hereby granted, provided that the above copyright + notice appear in all copies. The author and McGill University + make no representations about the suitability of this + software for any purpose. It is provided "as is" without + express or implied warranty. +---------------------------------------------------------------------------- */ + + +#include +#include +#include +#include +#include +#include +#include +#include "ecat_write.h" +#define ShortMax 32767 +#define BECQUEREL_PER_NCURIE 37 + +typedef enum MODALITY { + UNKNOWN_MODALITY, + PET, + MRI, + NUM_MODALITY +} modality; + +typedef enum DATA_UNIT { + UNKNOWN_DATA_UNIT, + LABEL, + NCIPERCC +} data_unit; + + +typedef struct MINC_INFO { + int dim; + long x_size; + long y_size; + long z_size; + long time_length; + double x_step; + double y_step; + double z_step; + int x_start_flag; + int y_start_flag; + int z_start_flag; + double x_start; + double y_start; + double z_start; + double *time_points; + double *time_widths; + data_unit dunit; + modality mod; + float isotope_halflife; +} minc_info; + + +/*Fuction declarations*/ +double decay_correction(double scan_time, double measure_time, double start_time, double half_life); +minc_info *get_minc_info(char *file); +void free_minc_info(minc_info *mi); +void get_patient_info(char *file, Main_header *mh); +void get_study_info(char *file, Main_header *mh); +void get_acquisition_info(char *file, Main_header *mh); +void get_ecat_acquisition_info(char *file, Main_header *mh); +void get_ecat_main_header_info(char *file, Main_header *mh); +void get_ecat_subheader_info(char *file, Image_subheader *sh); +Main_header * init_main_header(char *file, minc_info *mi); +Image_subheader * init_image_subheader(char *file, minc_info *mi, Main_header *mh); +void write_ecat_frame(MatrixFile *mf, Image_subheader *sh, short *ptr,int fr, float min, float max); +void usage_error(char *progname); + +int main ( int argc, char **argv) { + /* Variables for arguments */ + static int include_patient_info = TRUE; + static int include_study_info = TRUE; + static int include_acquisition_info = TRUE; + static int include_ecat_acquisition_info = TRUE; + static int include_ecat_main_info = TRUE; + static int include_ecat_subheader_info = TRUE; + static int do_decay_corr_factor = TRUE; + static int label_data = FALSE; + + /* Argument option table */ + static ArgvInfo argTable[] = { + {"-ignore_patient_variable", ARGV_CONSTANT, (char *) FALSE, (char *) &include_patient_info,"Ignore informations from the minc patient variable."}, + {"-ignore_study_variable", ARGV_CONSTANT, (char *) FALSE, (char *) &include_study_info,"Ignore informations from the minc study variable."}, + {"-ignore_acquisition_variable", ARGV_CONSTANT, (char *) FALSE, (char *) &include_acquisition_info,"Ignore informations from the minc acquisition variable."}, + {"-ignore_ecat_acquisition_variable", ARGV_CONSTANT, (char *) FALSE, (char *) &include_ecat_acquisition_info,"Ignore informations from the minc ecat_acquisition variable."}, + {"-ignore_ecat_main", ARGV_CONSTANT, (char *) FALSE, (char *) &include_ecat_main_info,"Ignore informations from the minc ecat-main variable."}, + {"-ignore_ecat_subheader_variable", ARGV_CONSTANT, (char *) FALSE, (char *) &include_ecat_subheader_info,"Ignore informations from the minc ecat-subhdr variable."}, + {"-no_decay_corr_fctr", ARGV_CONSTANT, (char *) FALSE, (char *) &do_decay_corr_factor,"Do not compute the decay correction factors"}, + {"-label", ARGV_CONSTANT, (char *) TRUE, (char *) &label_data,"Voxel values are treated as integers, scale and calibration factors are set to unity"}, + {NULL, ARGV_END, NULL, NULL, NULL} + }; + + /*other variables*/ + char *pname; /*name of the present command ->argv[0]*/ + char *minc_file_name = NULL; /*name of the input file (minc format)*/ + char *ecat_file_name = NULL; /*name of the output file (to be created into the ecat format*/ + minc_info *mi = NULL; + Volume volume; + short int *ecat_short_ptr = NULL; + float ecat_scale_factor = 1; + progress_struct progress; + STRING DimOrder[MAX_DIMENSIONS]; + MatrixFile *out_matrix_file = NULL; + Main_header *out_main_header = NULL; + Image_subheader *out_image_subheader = NULL; + int fr; + int incre = 0; + int x_ind, y_ind, z_ind; + double cal_factor; + + /* Check arguments */ + pname = argv[0]; + if (ParseArgv(&argc, argv, argTable, 0) || (argc != 3)) { + usage_error(pname); + } + + /* Get file names */ + minc_file_name = argv[1]; + ecat_file_name = argv[2]; + + /*getting informations about the dimensions from the minc file*/ + if((mi = get_minc_info(minc_file_name)) == NULL) { + fprintf(stderr, "Can not get informations from %s.\n", minc_file_name); + return (1); + } + + /*checking the number of dimension*/ + if((mi->dim != 3) && (mi->dim != 4)) { + fprintf(stderr, "input file has not a valid number of dimension (%d), should be either 3 (static file) or 4 (dynamic file)"); + return (-1); + } + + + /*************************************************************************************** + * filling the ecat header fields with the available informations * + ***************************************************************************************/ + /*MAIN HEADER: filling compulsory fields*/ + if((out_main_header = init_main_header(minc_file_name, mi)) == NULL) { + fprintf(stderr, "Can not initilize the ecat main header.\n"); + return (1); + } + /*additional informations if available and requested*/ + /*getting patient informations from the MIpatient variable*/ + + if(include_patient_info) + get_patient_info(minc_file_name, out_main_header); + /*getting study information from the MIstudy variable*/ + if(include_study_info) + get_study_info(minc_file_name, out_main_header); + /*getting acquisition information from the MIacquisition variable*/ + if(include_acquisition_info) + get_acquisition_info(minc_file_name, out_main_header); + /*getting acquisition information from the ecat_acquisition variable*/ + if(include_ecat_acquisition_info) + get_ecat_acquisition_info(minc_file_name, out_main_header); + /*getting ecat main header informations from the ecat-main variable*/ + if(include_ecat_main_info) + get_ecat_main_header_info(minc_file_name, out_main_header); + + /*SUBHEADER TEMPLATE: filling compulsory fields*/ + if((out_image_subheader = init_image_subheader(minc_file_name, mi, out_main_header)) == NULL) { + fprintf(stderr, "Can not initialize the ecat image subheader.\n"); + return (1); + } + /*getting additionnal informations from the ecat-subhdr variable*/ + if(include_ecat_subheader_info) + get_ecat_subheader_info(minc_file_name, out_image_subheader); + + + /*writting the main header*/ + if((out_matrix_file = matrix_create(ecat_file_name, out_main_header)) == NULL) { + fprintf(stderr, "cannot create %s \n", ecat_file_name); + return(-1); + } + + /*extracting the mincvolume in the order time, X, Y, Z*/ + if(mi->dim == 4) { + DimOrder[0] = MItime; + DimOrder[1] = MIzspace; + DimOrder[2] = MIyspace; + DimOrder[3] = MIxspace; + } else { + DimOrder[0] = MIzspace; + DimOrder[1] = MIyspace; + DimOrder[2] = MIxspace; + } + + + /*input the volume*/ + if(input_volume(minc_file_name, + mi->dim, + DimOrder, + NC_UNSPECIFIED, + FALSE, + 1, + 1, + TRUE, + &volume, + (minc_input_options *) NULL) != OK) return (1); + + /*initializing the progress report*/ + initialize_progress_report(&progress, + FALSE, + mi->time_length, + "Converting data"); + cal_factor = (mi->dunit == NCIPERCC)?BECQUEREL_PER_NCURIE:1; + cal_factor /= out_main_header->calibration_factor; + + /*getting the values*/ + for(fr = 0; fr < ((mi->time_length)?mi->time_length:1); fr++) { + float max_val = 0, min_val = 0; + int vx, vy, vz; + Real value; + + if((ecat_short_ptr = (short int *)calloc((mi->x_size) * (mi->y_size) * (mi->z_size), sizeof(short int))) == NULL) { + fprintf(stderr, "allocation error\n"); + return 1; + } + + for(vz = 0; vz < mi->z_size; vz++) + for(vy = 0; vy < mi->y_size; vy++) + for(vx = 0; vx < mi->x_size; vx++) { + value = ((mi->time_length)?get_volume_real_value(volume, fr, vz, vy, vx, 0):get_volume_real_value(volume, vz, vy, vx, 0,0)) * cal_factor; + if(value > max_val) max_val = value; + if(value < min_val) min_val = value; + } + if((mi->dunit != LABEL) && (!label_data)) + ecat_scale_factor = ((max_val > fabs(min_val))?max_val:fabs(min_val))/ShortMax; + + for(vz = 0; vz < mi->z_size; vz++) { + for(vy = 0; vy < mi->y_size; vy++) { + for(vx = 0; vx < mi->x_size; vx++) { + value = ((mi->time_length)?get_volume_real_value(volume, fr, vz, vy, vx, 0):get_volume_real_value(volume, vz, vy, vx, 0,0)); + x_ind = ((mi->x_step > 0)?mi->x_size - 1 - vx:vx); + y_ind = ((mi->y_step > 0)?mi->y_size - 1 - vy:vy); + z_ind = ((mi->z_step > 0)?mi->z_size - 1 - vz:vz); + ecat_short_ptr[IJK(z_ind, y_ind, x_ind, mi->y_size, mi->x_size)] = (short int)ROUND(value * cal_factor/(ecat_scale_factor)); + } + } + } + + /*filling image subheader*/ + if(mi->dim == 4) { + out_image_subheader->frame_start_time = (unsigned int)(mi->time_points[fr] * 1000); + out_image_subheader->frame_duration = (unsigned int)(mi->time_widths[fr] * 1000); + } + if((mi->dunit == NCIPERCC) && (do_decay_corr_factor)) + out_image_subheader->decay_corr_fctr = decay_correction(mi->time_points[fr],mi->time_widths[fr] , 0.0, mi->isotope_halflife); + out_image_subheader->scale_factor = ecat_scale_factor; + write_ecat_frame(out_matrix_file, out_image_subheader, ecat_short_ptr, fr, min_val, max_val); + free(ecat_short_ptr); + ecat_short_ptr = NULL; + update_progress_report(&progress, incre); + incre++; + } + free_minc_info(mi); + matrix_close(out_matrix_file); + free(out_main_header); + free(out_image_subheader); + + terminate_progress_report(&progress); + return (0); +} + + + +double decay_correction(double scan_time, double measure_time, double start_time, double half_life) { + double mean_life; + double measure_correction; + double decay; + + /* Check for negative half_life and calculate mean life */ + if (half_life <= 0.0) return 1.0; + mean_life = half_life/ log(2.0); + + /* Normalize scan time and measure_time */ + scan_time = (scan_time - start_time) / mean_life; + measure_time /= mean_life; + + /* Calculate correction for decay over measuring time (assuming a + constant activity). Check for possible rounding errors. */ + if ((measure_time*measure_time/2.0) < DBL_EPSILON) { + measure_correction = 1.0 - measure_time/2.0; + } + else { + measure_correction = (1.0 - exp(-measure_time)) / fabs(measure_time); + } + + /* Calculate decay */ + decay = exp(-scan_time) * measure_correction; + if (decay<=0.0) decay = DBL_MAX; + else decay = 1.0/decay; + + return decay; +} + +minc_info *get_minc_info(char *file) { + int minc_fd; /*minc file descriptor*/ + minc_info *mi = NULL; /*minc info structures*/ + int var_id; + int num_att; + char attname[NC_MAX_NAME+1]; + char buff[NC_MAX_NAME+1]; + int i; + + /* we first open the netcdf file*/ + minc_fd = miopen(file, NC_NOWRITE); + + if((mi = (minc_info *)calloc(1, sizeof(minc_info))) == NULL) return NULL; + + /*initlaization*/ + mi->time_length = 0; + mi->mod = UNKNOWN_MODALITY; + mi->dunit = UNKNOWN_DATA_UNIT; + mi->x_start_flag = 0; + mi->y_start_flag = 0; + mi->z_start_flag = 0; + + /* first we check out the number of dimension. we suppose that the time dimension is the 4th*/ + ncinquire(minc_fd, &(mi->dim), NULL, NULL, NULL); + + /*getting infos about the x, y and and z sizes and width*/ + ncdiminq(minc_fd, ncdimid(minc_fd, MIxspace), NULL, &(mi->x_size)); + ncdiminq(minc_fd, ncdimid(minc_fd, MIyspace), NULL, &(mi->y_size)); + ncdiminq(minc_fd, ncdimid(minc_fd, MIzspace), NULL, &(mi->z_size)); + + /*getting information about each dimensions*/ + /*MIxspace*/ + var_id = ncvarid(minc_fd, MIxspace); + ncvarinq(minc_fd,var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + ncattname(minc_fd,var_id,i,attname); + if(strcmp(MIstart, attname) == 0) { + mi->x_start_flag = 1; + miattget1(minc_fd,var_id, MIstart, NC_DOUBLE, &(mi->x_start)); + } + if(strcmp(MIstep, attname) == 0) + miattget1(minc_fd,var_id, MIstep, NC_DOUBLE, &(mi->x_step)); + } + /*MIyspace*/ + var_id = ncvarid(minc_fd, MIyspace); + ncvarinq(minc_fd,var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + ncattname(minc_fd,var_id,i,attname); + if(strcmp(MIstart, attname) == 0) { + mi->y_start_flag = 1; + miattget1(minc_fd,var_id, MIstart, NC_DOUBLE, &(mi->y_start)); + } + if(strcmp(MIstep, attname) == 0) + miattget1(minc_fd,var_id, MIstep, NC_DOUBLE, &(mi->y_step)); + } + + /*MIzspace*/ + var_id = ncvarid(minc_fd, MIzspace); + ncvarinq(minc_fd,var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + ncattname(minc_fd,var_id,i,attname); + if(strcmp(MIstart, attname) == 0) { + mi->z_start_flag = 1; + miattget1(minc_fd,var_id, MIstart, NC_DOUBLE, &(mi->z_start)); + } + if(strcmp(MIstep, attname) == 0) + miattget1(minc_fd,var_id, MIstep, NC_DOUBLE, &(mi->z_step)); + } + + /* if it is a dynamic file, then we ask for the number of frames, the frame durations and start time*/ + if(mi->dim == 4){ + long start_time_vector[1]; + long count_time_vector[1]; + + var_id = ncdimid(minc_fd, MItime); + ncdiminq(minc_fd, var_id, NULL, &(mi->time_length)); + if((mi->time_points = (double *) calloc(mi->time_length, sizeof(double))) == NULL) return NULL; + if((mi->time_widths = (double *) calloc(mi->time_length, sizeof(double))) == NULL) return NULL; + start_time_vector[0] = 0; + count_time_vector[0] = mi->time_length; + + mivarget(minc_fd, var_id, start_time_vector, count_time_vector, NC_DOUBLE, MI_SIGNED, mi->time_points); + mivarget(minc_fd, ncvarid(minc_fd, MItime_width), start_time_vector, count_time_vector, NC_DOUBLE, MI_SIGNED, mi->time_widths); + } + + /*defining the modality*/ + if(mivar_exists(minc_fd, MIstudy)) { + var_id = ncvarid(minc_fd, MIstudy); + ncvarinq(minc_fd, var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + ncattname(minc_fd,var_id,i,attname); + if(strcmp(MImodality, attname) == 0) { + miattgetstr(minc_fd, var_id, MImodality, NC_MAX_NAME, buff); + if(strcmp(buff, MI_PET) == 0) mi->mod = PET; + if(strcmp(buff, MI_MRI) == 0) mi->mod = MRI; + if(strcmp(buff, MI_LABEL) == 0) mi->dunit = LABEL; + } + } + } + if(mivar_exists(minc_fd, MIimagemin)) { + var_id = ncvarid(minc_fd, MIimagemin); + ncvarinq(minc_fd, var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + ncattname(minc_fd,var_id,i,attname); + if(strcmp(MIunits, attname) == 0) { + miattgetstr(minc_fd, var_id, MIunits, NC_MAX_NAME, buff); + if(strcmp(buff, "nCi/cc") == 0) mi->dunit = NCIPERCC; + } + } + } + + /*getting isotope halflife if available*/ + if(mivar_exists(minc_fd, MIacquisition)) { + var_id = ncvarid(minc_fd, MIacquisition); + ncvarinq(minc_fd, var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + ncattname(minc_fd,var_id,i,attname); + if(strcmp(MIradionuclide_halflife, attname) == 0) + miattget1(minc_fd, var_id, MIradionuclide_halflife, NC_FLOAT, &(mi->isotope_halflife)); + } + } + + miclose(minc_fd); + return mi; +} + +void free_minc_info(minc_info *mi) { + if(mi->time_points != NULL) free(mi->time_points); + if(mi->time_widths != NULL) free(mi->time_widths); +} +void get_patient_info(char *file, Main_header *mh) { + int minc_fd; /*minc file descriptor*/ + + /* we first open the netcdf file*/ + minc_fd = miopen(file, NC_NOWRITE); + + if(mivar_exists(minc_fd, MIpatient)) { + int var_id; + int num_att; + int i; + + var_id = ncvarid(minc_fd, MIpatient); + ncvarinq(minc_fd, var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + char attname[NC_MAX_NAME+1]; + char buffer_val[NC_MAX_NAME+1]; + + ncattname(minc_fd,var_id,i,attname); + if(strcmp(MIfull_name, attname) == 0) + miattgetstr(minc_fd, var_id, MIfull_name, 32, mh->patient_name); + if(strcmp(MIage, attname) == 0) + miattget1(minc_fd, var_id, MIage, NC_FLOAT, &(mh->patient_age)); + if(strcmp(MIweight, attname) == 0) + miattget1(minc_fd, var_id, MIweight, NC_FLOAT, &(mh->patient_weight)); + if(strcmp(MIsize, attname) == 0) + miattget1(minc_fd, var_id, MIsize, NC_FLOAT, &(mh->patient_height)); + if(strcmp(MIsex, attname) == 0) { + miattgetstr(minc_fd, var_id, MIsex, NC_MAX_NAME+1, buffer_val); + if(strcmp(buffer_val, MI_MALE) == 0) + mh->patient_sex[0] = 'M'; + if(strcmp(buffer_val, MI_FEMALE) == 0) + mh->patient_sex[0] = 'F'; + } + if(strcmp(MIidentification, attname) == 0) + miattgetstr(minc_fd, var_id, MIidentification, 16, mh->patient_id); + + /*I don't kown how to deal with the patient birth date yet + in minc, it is a string, in ecat a int ... + I guess this is not important as long as we have the age of the patient + Might be filled below if the ecat-header variable exists + */ + } + } + miclose(minc_fd); +} + +void get_study_info(char *file, Main_header *mh) { + int minc_fd; /*minc file descriptor*/ + + /* we first open the netcdf file*/ + minc_fd = miopen(file, NC_NOWRITE); + + if(mivar_exists(minc_fd, MIstudy)) { + int var_id; + int i; + struct tm scan_time; + int num_att; + + var_id = ncvarid(minc_fd, MIstudy); + ncvarinq(minc_fd, var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + char attname[NC_MAX_NAME+1]; + + ncattname(minc_fd,var_id,i,attname); + if(strcmp(MIstart_year, attname) == 0) { + miattget1(minc_fd, var_id, MIstart_year, NC_INT, &scan_time.tm_year); + scan_time.tm_year -= 1900; + } + if(strcmp(MIstart_month, attname) == 0) { + miattget1(minc_fd, var_id, MIstart_month, NC_INT, &scan_time.tm_mon); + scan_time.tm_mon--; + } + if(strcmp(MIstart_day, attname) == 0) + miattget1(minc_fd, var_id, MIstart_day, NC_INT, &scan_time.tm_mday); + + if(strcmp(MIstart_hour, attname) == 0) + miattget1(minc_fd, var_id, MIstart_hour, NC_INT, &scan_time.tm_hour); + + if(strcmp(MIstart_minute, attname) == 0) + miattget1(minc_fd, var_id, MIstart_minute, NC_INT, &scan_time.tm_min); + + if(strcmp(MIstart_seconds, attname) == 0) + miattget1(minc_fd, var_id, MIstart_seconds, NC_INT, &scan_time.tm_sec); + + if(strcmp(MIidentification, attname) == 0) + miattgetstr(minc_fd, var_id, MIidentification, 16, mh->patient_id); + + if(strcmp(MIstudy_id, attname) == 0) + miattgetstr(minc_fd, var_id, MIstudy_id, 12, mh->study_name); + + if(strcmp(MIreferring_physician, attname) == 0) + miattgetstr(minc_fd, var_id, MIreferring_physician, 32, mh->physician_name); + + if(strcmp(MIoperator, attname) == 0) + miattgetstr(minc_fd, var_id, MIoperator, 32, mh->operator_name); + } + mh->scan_start_time = mktime(&scan_time); + } + miclose(minc_fd); +} + +void get_acquisition_info(char *file, Main_header *mh) { + int minc_fd; /*minc file descriptor*/ + + /* we first open the netcdf file*/ + minc_fd = miopen(file, NC_NOWRITE); + + if(mivar_exists(minc_fd, MIacquisition)) { + int var_id; + int num_att; + int i; + struct tm injection_time; + + var_id = ncvarid(minc_fd, MIacquisition); + ncvarinq(minc_fd, var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + char attname[NC_MAX_NAME+1]; + char buffer_val[NC_MAX_NAME+1]; + + ncattname(minc_fd,var_id,i,attname); + + if(strcmp(MIradionuclide, attname) == 0) + miattgetstr(minc_fd, var_id, MIradionuclide, 8, mh->isotope_code); + + if(strcmp(MIradionuclide_halflife, attname) == 0) + miattget1(minc_fd, var_id, MIradionuclide_halflife, NC_FLOAT, &(mh->isotope_halflife)); + + if(strcmp(MItracer, attname) == 0) + miattgetstr(minc_fd, var_id, MItracer, 32, mh->radiopharmaceutical); + + if(strcmp(MIinjection_year, attname) == 0) { + miattget1(minc_fd, var_id, MIinjection_year, NC_INT, &injection_time.tm_year); + injection_time.tm_year -= 1900; + } + if(strcmp(MIinjection_month, attname) == 0) { + miattget1(minc_fd, var_id, MIinjection_month, NC_INT, &injection_time.tm_mon); + injection_time.tm_mon--; + } + if(strcmp(MIinjection_day, attname) == 0) + miattget1(minc_fd, var_id, MIinjection_day, NC_INT, &injection_time.tm_mday); + + if(strcmp(MIinjection_hour, attname) == 0) + miattget1(minc_fd, var_id, MIinjection_hour, NC_INT, &injection_time.tm_hour); + + if(strcmp(MIinjection_minute, attname) == 0) + miattget1(minc_fd, var_id, MIinjection_minute, NC_INT, &injection_time.tm_min); + + if(strcmp(MIinjection_seconds, attname) == 0) + miattget1(minc_fd, var_id, MIinjection_seconds, NC_INT, &injection_time.tm_sec); + + if(strcmp(MIinjection_dose, attname) == 0) + miattget1(minc_fd, var_id, MIinjection_dose, NC_FLOAT, &(mh->dosage)); + + if(strcmp(MIdose_units, attname) == 0) + miattgetstr(minc_fd, var_id, MIdose_units, 32, mh->data_units); + } + mh->dose_start_time = mktime(&injection_time); + } + miclose(minc_fd); +} + +void get_ecat_acquisition_info(char *file, Main_header *mh) { + int minc_fd; /*minc file descriptor*/ + + /* we first open the netcdf file*/ + minc_fd = miopen(file, NC_NOWRITE); + + if(mivar_exists(minc_fd, "ecat_acquisition")) { + int var_id; + int num_att; + int i; + char attname[NC_MAX_NAME+1]; + + var_id = ncvarid(minc_fd, "ecat_acquisition"); + ncvarinq(minc_fd, var_id, NULL, NULL, NULL, NULL, &num_att); + for(i = 0; i < num_att; i++) { + ncattname(minc_fd,var_id,i,attname); + if(strcmp("septa_retracted", attname) == 0) { + char buff[NC_MAX_NAME+1]; + + miattgetstr(minc_fd, var_id, "septa_retracted", NC_MAX_NAME+1, buff); + if(strcmp(buff, MI_TRUE) == 0) + mh->septa_state = 1; + if(strcmp(buff, MI_FALSE) == 0) + mh->septa_state = 0; + } + } + } + miclose(minc_fd); +} + +void get_ecat_main_header_info(char *file, Main_header *mh) { + int minc_fd; + + /* we first open the netcdf file*/ + minc_fd = miopen(file, NC_NOWRITE); + + if(mivar_exists(minc_fd, "ecat-main")) { + int var_id; + char buffer[NC_MAX_NAME+1]; + + + var_id = ncvarid(minc_fd, "ecat-main"); + + miattgetstr(minc_fd, var_id, "Magic_Number", 14, mh->magic_number); + + miattgetstr(minc_fd, var_id, "Original_Filename", 32, mh->original_file_name); + + miattgetstr(minc_fd, var_id, "System_Type", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->system_type)); + + miattgetstr(minc_fd, var_id, "Serial_Number", 10, mh->serial_number); + + miattgetstr(minc_fd, var_id, "Scan_Start_Time", sizeof(buffer), buffer); + sscanf(buffer, "%u", &(mh->scan_start_time)); + + miattgetstr(minc_fd, var_id, "Isotope_Name", 8, mh->isotope_code); + + miattgetstr(minc_fd, var_id, "Isotope_Halflife", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->isotope_halflife)); + + miattgetstr(minc_fd, var_id, "Radiopharmaceutical", 32, mh->radiopharmaceutical); + + miattgetstr(minc_fd, var_id, "Gantry_Tilt", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->gantry_tilt)); + + miattgetstr(minc_fd, var_id, "Gantry_Rotation", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->gantry_rotation)); + + miattgetstr(minc_fd, var_id, "Bed_Elevation", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->bed_elevation)); + + miattgetstr(minc_fd, var_id, "Intrinsic_Tilt", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->intrinsic_tilt)); + + miattgetstr(minc_fd, var_id, "Wobble_Speed", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->wobble_speed)); + + miattgetstr(minc_fd, var_id, "Transm_Source_Type", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->transm_source_type)); + + miattgetstr(minc_fd, var_id, "Distance_Scanned", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->distance_scanned)); + + miattgetstr(minc_fd, var_id, "Transaxial_Fov", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->transaxial_fov)); + + miattgetstr(minc_fd, var_id, "Angular_Compression", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->angular_compression)); + + miattgetstr(minc_fd, var_id, "Coin_Samp_Mode", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->coin_samp_mode)); + + miattgetstr(minc_fd, var_id, "Axial_Samp_Mode", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->axial_samp_mode)); + + miattgetstr(minc_fd, var_id, "Calibration_Factor", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->calibration_factor)); + + miattgetstr(minc_fd, var_id, "Study_Type", 12, mh->study_name); + + miattgetstr(minc_fd, var_id, "Patient_Id", 16, mh->patient_id); + + miattgetstr(minc_fd, var_id, "Patient_Name", 32, mh->patient_name); + + miattgetstr(minc_fd, var_id, "Patient_Sex", sizeof(buffer), buffer); + sscanf(buffer, "%c", mh->patient_sex); + + miattgetstr(minc_fd, var_id, "Patient_Dexterity", sizeof(buffer), buffer); + sscanf(buffer, "%c", mh->patient_dexterity); + + miattgetstr(minc_fd, var_id, "Patient_Age", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->patient_age)); + + miattgetstr(minc_fd, var_id, "Patient_Height", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->patient_height)); + + miattgetstr(minc_fd, var_id, "Patient_Weight", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->patient_weight)); + + miattgetstr(minc_fd, var_id, "Patient_Birth_Date", sizeof(buffer), buffer); + sscanf(buffer, "%d", &(mh->patient_birth_date)); + + miattgetstr(minc_fd, var_id, "Physician_Name", 32, mh->physician_name); + + miattgetstr(minc_fd, var_id, "Operator_Name", 32, mh->operator_name); + + miattgetstr(minc_fd, var_id, "Study_Description", 32, mh->study_description); + + /*bugg -> in minc lib, acquisition should be replaced by acquisition*/ + miattgetstr(minc_fd, var_id, "Acquision_Type", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->acquisition_type)); + + miattgetstr(minc_fd, var_id, "Patient_Orientation", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->patient_orientation)); + + miattgetstr(minc_fd, var_id, "Facility_Name", 20, mh->facility_name); + + miattgetstr(minc_fd, var_id, "Lwr_Sctr_Thres", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->lwr_sctr_thres)); + + miattgetstr(minc_fd, var_id, "Lwr_True_Thres", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->lwr_true_thres)); + + miattgetstr(minc_fd, var_id, "Upr_True_Thres", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->upr_true_thres)); + + miattgetstr(minc_fd, var_id, "User_Process_Code", 10, mh->user_process_code); + + miattgetstr(minc_fd, var_id, "Acquisition_Mode", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->acquisition_mode)); + + miattgetstr(minc_fd, var_id, "Bin_Size", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->bin_size)); + + miattgetstr(minc_fd, var_id, "Branching_Fraction", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->branching_fraction)); + + miattgetstr(minc_fd, var_id, "Dose_Start_Time", sizeof(buffer), buffer); + sscanf(buffer, "%u", &(mh->dose_start_time)); + + miattgetstr(minc_fd, var_id, "Dosage", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->dosage)); + + miattgetstr(minc_fd, var_id, "Well_Counter_Corr_Factor", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(mh->well_counter_factor)); + + miattgetstr(minc_fd, var_id, "Data_Units", 32, mh->data_units); + + miattgetstr(minc_fd, var_id, "Septa_State", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(mh->septa_state)); + } + miclose(minc_fd); +} + +void get_ecat_subheader_info(char *file, Image_subheader *sh) { + int minc_fd; /*minc file descriptor*/ + + /* we first open the netcdf file*/ + minc_fd = miopen(file, NC_NOWRITE); + + if(mivar_exists(minc_fd, "ecat-subhdr")) { + int var_id; + char buffer[NC_MAX_NAME]; + + var_id = ncvarid(minc_fd, "ecat-subhdr"); + + miattgetstr(minc_fd, var_id, "Recon_Zoom", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->recon_zoom)); + + miattgetstr(minc_fd, var_id, "Filter_Code", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(sh->filter_code)); + + miattgetstr(minc_fd, var_id, "X_Resolution", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->x_resolution)); + + miattgetstr(minc_fd, var_id, "Y_Resolution", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->y_resolution)); + + miattgetstr(minc_fd, var_id, "Z_Resolution", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->z_resolution)); + + miattgetstr(minc_fd, var_id, "X_Rotation_Angle", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->num_r_elements)); + + miattgetstr(minc_fd, var_id, "Y_Rotation_Angle", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->num_angles)); + + miattgetstr(minc_fd, var_id, "Z_Rotation_Angle", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->z_rotation_angle)); + + miattgetstr(minc_fd, var_id, "Decay_Corr_Fctr", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->decay_corr_fctr)); + + miattgetstr(minc_fd, var_id, "Corrections_Applied", sizeof(buffer), buffer); + sscanf(buffer, "%d", &(sh->processing_code)); + + miattgetstr(minc_fd, var_id, "Gate_Duration", sizeof(buffer), buffer); + sscanf(buffer, "%u", &(sh->gate_duration)); + + miattgetstr(minc_fd, var_id, "R_Wave_Offset", sizeof(buffer), buffer); + sscanf(buffer, "%d", &(sh->r_wave_offset)); + + miattgetstr(minc_fd, var_id, "Num_Accepted_Beats", sizeof(buffer), buffer); + sscanf(buffer, "%d", &(sh->num_accepted_beats)); + + miattgetstr(minc_fd, var_id, "Filter_Cutoff_Frequency", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->filter_cutoff_frequency)); + + miattgetstr(minc_fd, var_id, "Filter_Dc_Component", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->filter_resolution)); + + miattgetstr(minc_fd, var_id, "Filter_Ramp_Slope", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->filter_ramp_slope)); + + miattgetstr(minc_fd, var_id, "Filter_Order", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(sh->filter_order)); + + miattgetstr(minc_fd, var_id, "Annotation", 40, sh->annotation); + + miattgetstr(minc_fd, var_id, "Da_X_Rotation_Angle", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->mt_1_1)); + + miattgetstr(minc_fd, var_id, "Da_Y_Rotation_Angle", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->mt_1_2)); + + miattgetstr(minc_fd, var_id, "Da_Z_Rotation_Angle", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->mt_1_3)); + + miattgetstr(minc_fd, var_id, "Da_X_Translation", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->mt_2_1)); + + miattgetstr(minc_fd, var_id, "Da_Y_Translation", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->mt_2_2)); + + miattgetstr(minc_fd, var_id, "Da_Z_Translation", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->mt_2_3)); + + miattgetstr(minc_fd, var_id, "Da_X_Scale_Factor", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->mt_3_1)); + + miattgetstr(minc_fd, var_id, "Da_Y_Scale_Factor", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->mt_3_2)); + + miattgetstr(minc_fd, var_id, "Da_Z_Scale_Factor", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->mt_3_3)); + + miattgetstr(minc_fd, var_id, "Rfilter_Cutoff", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->rfilter_cutoff)); + + miattgetstr(minc_fd, var_id, "Rfilter_Resolution", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->rfilter_resolution)); + + miattgetstr(minc_fd, var_id, "Rfilter_Code", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(sh->rfilter_code)); + + miattgetstr(minc_fd, var_id, "Rfilter_Order", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(sh->rfilter_order)); + + miattgetstr(minc_fd, var_id, "Zfilter_Cutoff", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->zfilter_cutoff)); + + miattgetstr(minc_fd, var_id, "Zfilter_Resolution", sizeof(buffer), buffer); + sscanf(buffer, "%f", &(sh->zfilter_resolution)); + + miattgetstr(minc_fd, var_id, "Zfilter_Code", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(sh->zfilter_code)); + + miattgetstr(minc_fd, var_id, "Zfilter_Order", sizeof(buffer), buffer); + sscanf(buffer, "%hd", &(sh->zfilter_order)); + } + miclose(minc_fd); +} + +Main_header * init_main_header(char *file, minc_info *mi) { + Main_header *mh = NULL; + + if((mh = (Main_header *)calloc(1, sizeof(Main_header))) == NULL) return NULL; + + /**************************************************************************** + * first, we filled the compulsory fields * + ****************************************************************************/ + + /*software version 7*/ + mh->sw_version = 72; + + /*well ecat is a format for PET data, so MRI data in the minc format will have the PetVolume flag in ecat ....*/ + mh->file_type = PetVolume; + + mh->num_planes = mi->z_size; + + mh->num_frames = (mi->dim == 4)?mi->time_length:1; + + mh->num_gates = 1; /*not handled for now*/ + + mh->num_bed_pos = 0; /*not handles for now*/ + + /*inital bed position, I do the backward calculation of ecattominc - !! MIstart may not exist*/ + mh->init_bed_position = -mi->z_start/10.; + mh->init_bed_position = (mi->z_step < 0)?-mi->z_start/10.:-(mi->z_step * (mi->z_size - 1) + mi->z_start)/10.; + mh->bed_offset[0] = 0; + + /*for now, the distance_scanned is equal to the number of axial planes times the plane width*/ + mh->distance_scanned = fabs(mi->z_step) * (mi->z_size + 1)/10.0; /*cm*/ + + mh->plane_separation = fabs(mi->z_step)/10.; + + mh->calibration_factor = 1; /*forced*/ + + mh->calibration_units = 2; /*2 = processed, by default*/ + + mh->calibration_units_label = 0; /*by default*/ + + if(mi->dunit == NCIPERCC) { + mh->calibration_units = 1; + mh->calibration_units_label = 1; + } + mh->patient_sex[0] = 'U'; /*for now unknown*/ + + mh->patient_dexterity[0] = 'U'; /*unknown for now*/; + + return mh; +} + +Image_subheader * init_image_subheader(char *file, minc_info *mi, Main_header *mh) { + Image_subheader *sh = NULL; + + if((sh = (Image_subheader *)calloc(1, sizeof(Image_subheader))) == NULL) return NULL; + + sh->data_type = SunShort; /*forced*/ + sh->num_dimensions = 3; /*forced*/ + + sh->x_dimension = mi->x_size; + sh->y_dimension = mi->y_size; + sh->z_dimension = mi->z_size; + sh->x_pixel_size = fabs(mi->x_step)/10.0; /*in cm*/ + sh->y_pixel_size = fabs(mi->y_step)/10.0; + sh->z_pixel_size = fabs(mi->z_step)/10.0; + sh->x_offset = 0; + if(mi->x_start_flag) + sh->x_offset = -(mi->x_start + (mi->x_size - 1.0)/2.0 * mi->x_step)/10.; + sh->y_offset = 0; + if(mi->y_start_flag) + sh->y_offset = (mi->y_start + (mi->y_size - 1.0)/2.0 * mi->y_step)/10.; + sh->z_offset = 0; + return sh; +} +void write_ecat_frame(MatrixFile *mf, Image_subheader *sh, short *ptr,int fr, float min, float max) { + MatrixData md = {0}; + + md.matfile = mf; + md.mat_type = mf->mhptr->file_type; + md.data_type = sh->data_type; + md.shptr = (caddr_t)sh; + md.data_ptr = (caddr_t)ptr; + md.xdim = sh->x_dimension; + md.ydim = sh->y_dimension; + md.zdim = sh->z_dimension; + + md.scale_factor = sh->scale_factor; + md.pixel_size = sh->x_pixel_size; + md.y_size = sh->y_pixel_size; + md.z_size = sh->z_pixel_size; + md.data_min = min; + md.data_max = max; + md.data_size = (sh->x_dimension) * (sh->y_dimension) * (sh->z_dimension) * sizeof(short int); + matrix_write(mf,mat_numcod(fr + 1, 1, 1, 0, 0) , &md); +} + + +void usage_error(char *progname) +{ + (void) fprintf(stderr, + "\nUsage: %s [] \n", progname); + (void) fprintf(stderr, + " %s [-help]\n\n", progname); + + exit(EXIT_FAILURE); +} + + +