/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * filename: m-nifti.c * * * * UTIL C-source: Medical Image Conversion Utility * * * * purpose : read and write NIFTI files * * * * project : (X)MedCon by Erik Nolf * * * * Functions : MdcCheckNIFTI() - Check NIFTI format * * MdcReadNIFTI() - Read NIFTI file * * MdcWriteNIFTI() - Write NIFTI file * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /* $Id: m-nifti.c,v 1.20 2007/09/11 13:24:36 enlf Exp $ */ /* Copyright (C) 1997-2007 by Erik Nolf This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Place - Suite 330, Boston, MA 02111-1307, USA. */ /**************************************************************************** H E A D E R S ****************************************************************************/ #include "m-depend.h" #include /* #include #ifdef HAVE_STDLIB_H #include #endif #ifdef HAVE_STRING_H #include #endif #ifdef HAVE_STRINGS_H #ifndef _WIN32 #include #endif #endif #ifdef HAVE_UNISTD_H #include #endif */ #include "medcon.h" /**************************************************************************** D E F I N E S *****************************************************************************/ #define MDC_NIFTI_WRITE_QFORM 0 /* 0/1 write qform orientation & location */ /**************************************************************************** F U N C T I O N S *****************************************************************************/ int MdcCheckNIFTI(FILEINFO *fi) { int ret, FORMAT=MDC_FRMT_NONE; MdcMergePath(fi->ipath,fi->idir,fi->ifname); nifti_set_debug_level(0); ret = is_nifti_file(fi->ipath); nifti_set_debug_level(1); MdcSplitPath(fi->ipath,fi->idir,fi->ifname); switch (ret) { #if MDC_INCLUDE_ANLZ case 0: FORMAT = MDC_FRMT_NONE; /* check later as Analyze */ break; #else case 0: FORMAT = MDC_FRMT_NIFTI; /* use NIFTI reader */ break; #endif case 1: FORMAT = MDC_FRMT_NIFTI; /* NIFTI one file */ break; case 2: FORMAT = MDC_FRMT_NIFTI; /* NIFTI two files */ break; default: FORMAT = MDC_FRMT_NONE; /* unknown */ } return(FORMAT); } const char *MdcReadNIFTI(FILEINFO *fi) { nifti_1_header *nhdr; nifti_image *nim; Uint8 *pbuf=NULL; char *fname; int swap, dim[8]; float xyz_fctr, t_fctr; Uint32 i, number, f, bytes; IMG_DATA *id=NULL; DYNAMIC_DATA *dd=NULL; if (MDC_FILE_STDIN == MDC_YES) return("NIFTI File input from stdin unsupported"); if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_BEGIN,0.,"Reading NIFTI:"); if (MDC_VERBOSE) MdcPrntMesg("NIFTI Reading <%s> ...",fi->ifname); MDC_FILE_ENDIAN = MDC_HOST_ENDIAN; MdcMergePath(fi->ipath,fi->idir,fi->ifname); if (MDC_INFO) { nhdr = nifti_read_header(fi->ipath,&swap,MDC_YES); if (nhdr == NULL) { MdcSplitPath(fi->ipath,fi->idir,fi->ifname); return("NIFTI Failure reading header"); } disp_nifti_1_header("NIFTI", nhdr); MdcFree(nhdr); } fname = malloc(strlen(fi->ipath)+4); if (fname == NULL) { MdcSplitPath(fi->ipath,fi->idir,fi->ifname); return("NIFTI Failure to malloc filename path"); } strcpy(fname,fi->ipath); MdcSplitPath(fi->ipath,fi->idir,fi->ifname); if (fi->compression == MDC_GZIP) strcat(fname,".gz"); /* orig filename */ nim = nifti_image_read(fname, MDC_YES); MdcFree(fname); if (nim == NULL) { return("NIFTI Failure reading image"); } /* fill in FILEINFO */ fi->reconstructed = MDC_YES; fi->acquisition_type = MDC_ACQUISITION_TOMO; fi->endian = MDC_FILE_ENDIAN; MdcStringCopy(fi->study_descr,nim->descrip,80); if (MDC_ECHO_ALIAS == MDC_YES) { nifti_image_free(nim); MdcEchoAliasName(fi); return(NULL); } /* copy dim parameters */ for (i=0; i<8; i++) { fi->dim[i] = nim->dim[i]; fi->pixdim[i] = nim->pixdim[i]; if ((i > 0) && (fi->pixdim[i])) fi->pixdim[0] = i-1; } /* catch special case for single image */ if (fi->dim[0] == 2) fi->dim[0]=3; /* get unit rescale */ switch (nim->xyz_units) { case NIFTI_UNITS_METER : xyz_fctr = 1000.; break; case NIFTI_UNITS_MICRON: xyz_fctr = 1./1000.; break; case NIFTI_UNITS_MM : default : xyz_fctr = 1.; } switch (nim->time_units) { case NIFTI_UNITS_SEC : t_fctr = 1000.; break; case NIFTI_UNITS_USEC : t_fctr = 1./1000.; break; case NIFTI_UNITS_MSEC : default : t_fctr = 1.; } /* scale to internal units */ fi->pixdim[1] *= xyz_fctr; /* mm */ fi->pixdim[2] *= xyz_fctr; /* mm */ fi->pixdim[3] *= xyz_fctr; /* mm */ if (fi->dim[4] > 1) fi->pixdim[4] *= t_fctr; /* ms */ fi->mwidth = (Uint32) nim->nx; fi->mheight = (Uint32) nim->ny; for ( number=1, i=3; i<=nim->dim[0]; i++) number*=nim->dim[i]; if (number == 0) { nifti_image_free(nim); return("NIFTI No valid images specified"); } switch (nim->datatype) { case NIFTI_TYPE_UINT8 : fi->type=BIT8_U; fi->bits=8; break; case NIFTI_TYPE_INT16 : fi->type=BIT16_S; fi->bits=16; break; case NIFTI_TYPE_INT32 : fi->type=BIT32_S; fi->bits=32; break; case NIFTI_TYPE_FLOAT32 : fi->type=FLT32; fi->bits=32; break; case NIFTI_TYPE_FLOAT64 : fi->type=FLT64; fi->bits=64; break; case NIFTI_TYPE_RGB24 : fi->type=COLRGB; fi->bits=24; break; case NIFTI_TYPE_INT8 : fi->type=BIT8_S; fi->bits=8; break; case NIFTI_TYPE_UINT16 : fi->type=BIT16_U; fi->bits=16; break; case NIFTI_TYPE_UINT32 : fi->type=BIT32_U; fi->bits=32; break; case NIFTI_TYPE_INT64 : fi->type=BIT64_S; fi->bits=64; break; case NIFTI_TYPE_UINT64 : fi->type=BIT64_U; fi->bits=64; break; case NIFTI_TYPE_COMPLEX64 : case NIFTI_TYPE_FLOAT128 : case NIFTI_TYPE_COMPLEX128 : case NIFTI_TYPE_COMPLEX256 : default : nifti_image_free(nim); return("NIFTI Unsupported data type"); } /* read image data */ if (nifti_image_load(nim) < 0) { nifti_image_free(nim); return("NIFTI Failure loading data"); } /* get IMG_DATA structs */ if (!MdcGetStructID(fi,number)) { nifti_image_free(nim); return("NIFTI Bad malloc IMG_DATA structs"); } /* fill in IMG_DATA structs */ bytes = fi->mwidth * fi->mheight * MdcType2Bytes(fi->type); pbuf = (Uint8 *)nim->data; i=0; dim[0]=7; for (dim[7]=0; dim[7] < nim->nw; dim[7]++ ) /* nw */ for (dim[6]=0; dim[6] < nim->nv; dim[6]++ ) /* nv */ for (dim[5]=0; dim[5] < nim->nu; dim[5]++ ) /* nu */ for (dim[4]=0; dim[4] < nim->nt; dim[4]++ ) /* nt */ for (dim[3]=0; dim[3] < nim->nz; dim[3]++, i++, pbuf+=bytes) { /* nz */ if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_INCR,1./(float)fi->number,NULL); if (i == fi->number) { nifti_image_free(nim); return("NIFTI Internal ERRRO"); } id = &fi->image[i]; id->width = fi->mwidth; id->height= fi->mheight; id->bits = fi->bits; id->type = fi->type; id->quant_scale = nim->scl_slope; if (id->quant_scale == 0.) id->quant_scale = 1.; id->intercept = nim->scl_inter; id->pixel_xsize = fi->pixdim[1]; id->pixel_ysize = fi->pixdim[2]; id->slice_width = fi->pixdim[3]; id->slice_spacing = id->slice_width; if ( (id->buf=MdcGetImgBuffer(bytes)) == NULL) { nifti_image_free(nim); return("NIFTI Bad malloc image buffer"); } memcpy(id->buf,pbuf,bytes); } /* check some final FILEINFO entries */ if (fi->dim[4] > 1) { fi->acquisition_type = MDC_ACQUISITION_DYNAMIC; /* fill in dynamic data struct */ if (!MdcGetStructDD(fi,(unsigned)fi->dim[4])) { nifti_image_free(nim); return("NIFTI Couldn't malloc DYNAMIC_DATA structs"); } for (f=0; f < fi->dynnr; f++) { dd = &fi->dyndata[f]; dd->nr_of_slices = fi->dim[3]; dd->time_frame_delay = nim->toffset * t_fctr; dd->time_frame_duration = fi->pixdim[4]; dd->time_frame_start = f * dd->time_frame_duration + dd->time_frame_delay; } } nifti_image_free(nim); return(NULL); } const char *MdcWriteNIFTI(FILEINFO *fi) { struct nifti_1_header nhdr; nifti_image *nim; znzFile fp=NULL; char *bname; int i, n, ret, nvox, FREE; IMG_DATA *id; Int8 saved_norm_over_frames=MDC_NORM_OVER_FRAMES; Uint8 *buf=NULL, *maxbuf, *rgbbuf, grval; Uint32 size; Int16 type; if (XMDC_GUI == MDC_NO) { MdcDefaultName(fi,MDC_FRMT_NIFTI,fi->ofname,fi->ifname); } if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_BEGIN,0.,"Writing NIFTI:"); if (MDC_VERBOSE) MdcPrntMesg("NIFTI Writing <%s> ...",fi->ofname); if (MDC_FILE_STDOUT == MDC_YES) { return("NIFTI Writing to stdout currently unsupported"); } /* file endian */ if (MDC_WRITE_ENDIAN != MDC_HOST_ENDIAN) return("NIFTI Writing in different endianess yet unsupported"); /* get nifti_image struct */ nim = nifti_simple_init_nim(); if (nim == NULL) { return("NIFTI Couldn't init nifti_image struct"); } /* fill in header */ /* dimensions */ nim->ndim = fi->dim[0]; nim->nx = fi->dim[1]; nim->ny = fi->dim[2]; nim->nz = fi->dim[3]; nim->nt = fi->dim[4]; nim->nu = fi->dim[5]; nim->nv = fi->dim[6]; nim->nw = fi->dim[7]; for (i=0; i<8; i++) nim->dim[i] = fi->dim[i]; for (i=1,nvox=1; i<8; i++) nvox *= fi->dim[i]; nim->dx = fi->pixdim[1]; nim->dy = fi->pixdim[2]; nim->dz = fi->pixdim[3]; nim->dt = fi->pixdim[4]; nim->du = fi->pixdim[5]; nim->dv = fi->pixdim[6]; nim->dw = fi->pixdim[7]; for (i=0; i<8; i++) nim->pixdim[i] = fi->pixdim[i]; nim->xyz_units = NIFTI_UNITS_MM; nim->time_units = NIFTI_UNITS_MSEC; #if MDC_NIFTI_WRITE_QFORM /* orientation and location */ id = &fi->image[0]; nim->qform_code = NIFTI_XFORM_SCANNER_ANAT; nim->qoffset_x = - id->image_pos_pat[0]; nim->qoffset_y = - id->image_pos_pat[1]; nim->qoffset_z = + id->image_pos_pat[2]; nim->qto_xyz = nifti_make_orthog_mat44(- id->image_pos_pat[0], - id->image_pos_pat[1], + id->image_pos_pat[2], - id->image_pos_pat[3], - id->image_pos_pat[4], + id->image_pos_pat[5], 0,0,0); nifti_mat44_to_quatern(nim->qto_xyz, &nim->quatern_b, &nim->quatern_c, &nim->quatern_c, NULL,NULL,NULL,NULL,NULL,NULL,&nim->qfac); #endif /* single file output */ nim->nifti_type = 1; sprintf(nim->descrip,"%.35s",fi->study_descr); bname = nifti_makebasename(fi->opath); if (MDC_FILE_OVERWRITE == MDC_YES) { ret = nifti_set_filenames(nim,bname,0,1); }else{ ret = nifti_set_filenames(nim,bname,1,1); } free(bname); if (ret < 0) { nifti_image_free(nim); return("NIFTI Filename creation failed"); } /* output pixel type */ if (fi->map == MDC_MAP_PRESENT) { /* colored */ nim->datatype = NIFTI_TYPE_RGB24; nim->nbyper = 3; }else{ /* grayscale */ if (MDC_FORCE_INT != MDC_NO) { switch (MDC_FORCE_INT) { case BIT8_U : nim->datatype = NIFTI_TYPE_UINT8; nim->nbyper = 1; break; case BIT16_S: default : nim->datatype = NIFTI_TYPE_INT16; nim->nbyper = 2; } }else if (fi->diff_type) { nim->datatype = NIFTI_TYPE_INT16; nim->nbyper = 2; }else if (fi->diff_scale) { nim->datatype = NIFTI_TYPE_FLOAT32; nim->nbyper = 4; }else{ nim->nbyper = MdcType2Bytes(fi->type); switch ( fi->type ) { case BIT8_S : nim->datatype = NIFTI_TYPE_INT8; break; case BIT8_U : nim->datatype = NIFTI_TYPE_UINT8; break; case BIT16_S: nim->datatype = NIFTI_TYPE_INT16; break; case BIT16_U: nim->datatype = NIFTI_TYPE_UINT16; break; case BIT32_S: nim->datatype = NIFTI_TYPE_INT32; break; case BIT32_U: nim->datatype = NIFTI_TYPE_UINT32; break; case BIT64_S: nim->datatype = NIFTI_TYPE_INT64; break; case BIT64_U: nim->datatype = NIFTI_TYPE_UINT64; break; case FLT32 : nim->datatype = NIFTI_TYPE_FLOAT32; break; case FLT64 : nim->datatype = NIFTI_TYPE_FLOAT64; break; case ASCII: case BIT1 : default : nifti_image_free(nim); return("NIFTI Unsupported datatype"); } } } /* initialize output file */ fp = nifti_image_write_hdr_img(nim,2,"wb"); /* keep lowlevel header */ nhdr = nifti_convert_nim2nhdr(nim); /* free nifti image struct */ nifti_image_free(nim); if (fp == NULL) return("NIFTI Writing header data failed"); /* rescale over all images for */ /* a single slope/intercept */ MDC_NORM_OVER_FRAMES = MDC_NO; /* write (addapted) image data */ for (i=0; inumber; i++) { if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_INCR,1./(float)fi->number,NULL); id = &fi->image[i]; buf = id->buf; FREE = MDC_NO; type = id->type; if (fi->map != MDC_MAP_PRESENT) { /* grayscale */ if (MDC_FORCE_INT != MDC_NO) { switch (MDC_FORCE_INT) { case BIT8_U : buf = MdcGetImgBIT8_U(fi,i); type = BIT8_U; FREE=MDC_YES; break; case BIT16_S: buf = MdcGetImgBIT16_S(fi,i); type = BIT16_S; FREE=MDC_YES; break; default : buf = MdcGetImgBIT16_S(fi,i); type = BIT16_S; FREE=MDC_YES; } }else if (fi->diff_type) { switch (id->type) { case BIT16_S: buf = id->buf; type = BIT16_S; FREE=MDC_NO; break; default : buf = MdcGetImgBIT16_S(fi,i); type = BIT16_S; FREE=MDC_YES; } }else if (fi->diff_scale) { switch (id->type) { case FLT32 : buf = id->buf; type = FLT32; FREE=MDC_NO; break; default : buf = MdcGetImgFLT32(fi,i); type = FLT32; FREE=MDC_YES; break; } }else{ /* all (or most) types supported */ buf = id->buf; FREE=MDC_NO; type = id->type; } } if (buf == NULL) { znzclose(fp); return("NIFTI Bad malloc image buffer"); } if (fi->diff_size) { maxbuf = MdcGetResizedImage(fi, buf, type, i); if (FREE) MdcFree(buf); if (maxbuf == NULL) { znzclose(fp); return("NIFTI Bad malloc maxbuf"); } FREE=MDC_YES; }else{ maxbuf = buf; } size = fi->mwidth * fi->mheight * MdcType2Bytes(type); if (fi->map == MDC_MAP_PRESENT) { if (type == COLRGB) { /* true color */ if (nifti_write_buffer(fp,(void *)maxbuf,size) != size) { if (FREE) MdcFree(maxbuf); znzclose(fp); return("NIFTI Bad write RGB buffer"); } }else{ /* indexed */ rgbbuf = malloc(size * 3); if (rgbbuf == NULL) { if (FREE) MdcFree(maxbuf); znzclose(fp); return("NIFTI Bad mallox indexed buffer"); } /* make true color */ for (n=0; n < size; n += MdcType2Bytes(type)) { grval = (Uint8)MdcGetDoublePixel((Uint8 *)&maxbuf[n],type); rgbbuf[n*3 + 0] = fi->palette[grval * 3 + 0]; /* red */ rgbbuf[n*3 + 1] = fi->palette[grval * 3 + 1]; /* green */ rgbbuf[n*3 + 2] = fi->palette[grval * 3 + 2]; /* blue */ } if (FREE) MdcFree(maxbuf); maxbuf = rgbbuf; FREE = MDC_YES; size *= 3; /* RGB triplets */ if (nifti_write_buffer(fp,(void *)maxbuf,size) != size) { if (FREE) MdcFree(maxbuf); znzclose(fp); return("NIFTI Writing indexed buffer failed"); } } }else{ /* grayscale */ if (nifti_write_buffer(fp,(void *)maxbuf,size) != size) { if (FREE) MdcFree(maxbuf); znzclose(fp); return("NIFTI Bad write image buffer"); } } if (FREE) MdcFree(maxbuf); } /* update (rescaled) slope/intercept */ if (fi->image[0].rescaled == MDC_YES) { nhdr.scl_slope = fi->image[0].rescaled_slope; nhdr.scl_inter = fi->image[0].rescaled_intercept; }else{ nhdr.scl_slope = fi->image[0].rescale_slope; nhdr.scl_inter = fi->image[0].rescale_intercept; } /* rewrite updated header */ znzseek(fp, 0L, SEEK_SET); if (znzwrite(&nhdr,1,sizeof(nhdr),fp) != sizeof(nhdr)) { znzclose(fp); return("NIFTI Failure to update header"); } /* close file */ znzclose(fp); /* restore original value */ MDC_NORM_OVER_FRAMES = saved_norm_over_frames; /* finish */ return(NULL); }