/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * filename: m-inw.c * * * * UTIL C-source: Medical Image Conversion Utility * * * * purpose : read and write INW 1.0 files * * * * project : (X)MedCon by Erik Nolf * * * * Functions : MdcCheckINW() - Check INW format * * MdcReadINW() - Read INW file * * MdcWriteHeadStart() - Write Head_start * * MdcWriteHeadGen() - Write Head_gen * * MdcSkipHeadSpecs() - Skip Head_specs in file * * MdcWriteHeadSpecs() - Write Head_specs * * MdcWriteINW() - Write INW file * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ /* $Id: m-inw.c,v 1.36 2007/05/21 20:16:13 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 #ifdef HAVE_STDLIB_H #include #endif #ifdef HAVE_STRING_H #include #endif #ifdef HAVE_STRINGS_H #ifndef _WIN32 #include #endif #endif #include "medcon.h" /**************************************************************************** F U N C T I O N S ****************************************************************************/ int MdcCheckINW(FILEINFO *fi) { MDC_INW_HEAD_START hs; MDC_FILE_ENDIAN = MDC_LITTLE_ENDIAN; /* always */ if (fread((char *)&hs,1,MDC_INW_HEAD_START_SIZE,fi->ifp) != MDC_INW_HEAD_START_SIZE) return(MDC_BAD_READ); MdcSWAP(hs.mark); if (hs.mark != MDC_INW_SIG) return(MDC_FRMT_NONE); return(MDC_FRMT_INW); } char *MdcReadINW(FILEINFO *fi) { FILE *fp = fi->ifp; MDC_INW_HEAD_START hs; MDC_INW_HEAD_GEN hg; MDC_INW_HEAD_SPEC *hsp; IMG_DATA *id; Uint32 i, bytes, number; char *err=NULL; if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_BEGIN,0.,"Reading INW:"); if (MDC_VERBOSE) MdcPrntMesg("INW Reading <%s> ...",fi->ifname); if (MDC_ECHO_ALIAS == MDC_YES) { MdcEchoAliasName(fi); return(NULL); /* useless here */ } memset(&hs,0,MDC_INW_HEAD_START_SIZE); memset(&hg,0,MDC_INW_HEAD_GEN_SIZE); if (fread((char *)&hs,1,MDC_INW_HEAD_START_SIZE,fp) != MDC_INW_HEAD_START_SIZE) return("INW Bad read HeadStart struct"); /* put some defaults we use */ fi->reconstructed = MDC_YES; fi->acquisition_type = MDC_ACQUISITION_TOMO; fi->endian=MDC_FILE_ENDIAN=MDC_LITTLE_ENDIAN; MdcSWAP(hs.mark); MdcSWAP(hs.version); MdcSWAP(hs.size_header); MdcSWAP(hs.size_start); MdcSWAP(hs.size_gen); MdcSWAP(hs.size_spec); if (fread((char *)&hg,1,MDC_INW_HEAD_GEN_SIZE,fp) != MDC_INW_HEAD_GEN_SIZE) return("INW Bad read HeadGen struct"); MdcSWAP(hg.no); MdcSWAP(hg.sizeX); MdcSWAP(hg.sizeY); MdcSWAP(hg.pixel_type); MdcSWAP(hg.init_trans); MdcSWAP(hg.dummy1); MdcSWAP(hg.time); MdcSWAP(hg.scanner); MdcMakeIEEEfl(hg.max); MdcMakeIEEEfl(hg.min); MdcMakeIEEEfl(hg.decay_cst); MdcMakeIEEEfl(hg.pixel_size); if (MDC_INFO) { MdcPrntScrn("\nHEAD START (%d bytes)\n",MDC_INW_HEAD_START_SIZE); MdcPrintLine('-',MDC_HALF_LENGTH); MdcPrntScrn("mark : 0x%x\n",hs.mark); MdcPrntScrn("version : %hd.%hd\n",hs.version/256 ,hs.version%256); MdcPrntScrn("size_header : %hd\n",hs.size_header); MdcPrntScrn("size_start : %hd\n",hs.size_start); MdcPrntScrn("size_gen : %hd\n",hs.size_gen); MdcPrntScrn("size_spec : %hd\n",hs.size_spec); MdcPrntScrn("reserved : %.10s",hs.reserved); MdcPrntScrn("\n"); MdcPrntScrn("\nHEAD GEN (%d bytes)\n",MDC_INW_HEAD_GEN_SIZE); MdcPrintLine('-',MDC_HALF_LENGTH); MdcPrntScrn("number planes : %hd\n",hg.no); MdcPrntScrn("number columns: %hd\n",hg.sizeX); MdcPrntScrn("number rows : %hd\n",hg.sizeY); MdcPrntScrn("pixel type : %hd\n",hg.pixel_type); MdcPrntScrn("init transl : %hd [mm]\n",hg.init_trans); MdcPrntScrn("dummy1 : %hd\n",hg.dummy1); MdcPrntScrn("day : %.12s\n",hg.day); MdcPrntScrn("time : %d [sec]\n",hg.time); MdcPrntScrn("decay constant: %+e\n",hg.decay_cst); MdcPrntScrn("pixel size : %+e [mm]\n",hg.pixel_size); MdcPrntScrn("scaled max : %+e\n",hg.max); MdcPrntScrn("scaled min : %+e\n",hg.min); MdcPrntScrn("scanner type : %hd ",hg.scanner); switch (hg.scanner) { case EcatII: MdcPrntScrn("(EcatII)"); break; case EcatIV: MdcPrntScrn("(EcatIV)"); break; default : MdcPrntScrn("(Unknown)"); } MdcPrntScrn("\n"); MdcPrntScrn("recon type : %c ",hg.reconstruction); switch (hg.reconstruction) { case reconFBP : MdcPrntScrn("(Filtered Backprojection)"); break; case reconMaxLikFV: MdcPrntScrn("(Maximum likelihood F. Vermeulen)"); break; case reconMaxLik : MdcPrntScrn("(Maximum likelihood T. De Backer)"); break; case reconMaxPos : MdcPrntScrn("(Maximum a Posteriori)"); break; default : MdcPrntScrn("(Unknown)"); } MdcPrntScrn("\n"); MdcPrntScrn("recon version : %d\n",(int)hg.recon_version); MdcPrntScrn("reserved : %.24s\n",hg.reserved); } /* check some supported things */ if (hg.pixel_type != 2) return("INW Unsupported pixel type"); if (hs.version != (Int16)(MDC_INW_VERS_HIGH*256 + MDC_INW_VERS_LOW)) return("INW Unsupported version"); /* fill in the FILEINFO struct */ number = hg.no; if (number == 0 ) return("INW No valid images specified"); fi->mwidth = hg.sizeX; fi->mheight = hg.sizeY; fi->bits = hg.pixel_type * 8; switch (hg.pixel_type) { case 1: fi->type = BIT8_U; break; case 2: fi->type = BIT16_S; break; /* only this is supported */ case 4: fi->type = BIT32_S; break; #ifdef HAVE_8BYTE_INT case 8: fi->type = BIT64_S; break; #endif } fi->dim[0] = 3; fi->dim[1] = hg.sizeX; fi->dim[2] = hg.sizeY; fi->dim[3] = hg.no; fi->pixdim[0] = 3.; fi->pixdim[1] = hg.pixel_size; fi->pixdim[2] = hg.pixel_size; switch (hg.reconstruction) { case reconFBP : strcpy(fi->recon_method ,"Filtered Backprojection"); break; case reconMaxLikFV: strcpy(fi->recon_method ,"Maximum likelihood F. Vermeulen"); break; case reconMaxLik : strcpy(fi->recon_method ,"Maximum likelihood T. De Backer"); break; case reconMaxPos : strcpy(fi->recon_method ,"Maximum a Posteriori"); break; } /* now we read the spec/image headers */ if ((hsp = (MDC_INW_HEAD_SPEC *)malloc(MDC_INW_HEAD_SPEC_SIZE*number)) == NULL) return("INW Bad malloc HeadSpec structs"); memset(hsp,0,MDC_INW_HEAD_SPEC_SIZE*number); if (fread((char *)hsp,MDC_INW_HEAD_SPEC_SIZE,number,fi->ifp) != number) { MdcFree(hsp); return("INW Bad read HeadSpec structs"); } if (!MdcGetStructID(fi,number)) return("INW Bad malloc IMG_DATA structs"); for (i=0; inumber; i++) { if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_INCR,1./(float)fi->number,NULL); id = &fi->image[i]; MdcSWAP(hsp[i].time); MdcSWAP(hsp[i].trans); MdcSWAP(hsp[i].max); MdcSWAP(hsp[i].min); MdcMakeIEEEfl(hsp[i].cal_cst); if (MDC_INFO) { MdcPrntScrn("\n"); MdcPrntScrn("\nHEAD SPEC #%.3d (%d bytes)\n",i+1 ,MDC_INW_HEAD_SPEC_SIZE); MdcPrintLine('-',MDC_HALF_LENGTH); MdcPrntScrn("relative time : %d\n",hsp[i].time); MdcPrntScrn("calibr const : %+e\n",hsp[i].cal_cst); MdcPrntScrn("pixel max : %d\n",hsp[i].max); MdcPrntScrn("pixel min : %d\n",hsp[i].min); MdcPrntScrn("relat transl : %hd [mm]\n",hsp[i].trans); MdcPrntScrn("reserved : %.6s\n",hsp[i].reserved); } if (hsp[i].trans < 0.0 ) hsp[i].trans = - hsp[i].trans; /* positive */ /* fill in the IMG_DATA struct */ id->width = fi->mwidth; id->height = fi->mheight; id->bits = fi->bits; id->type = fi->type; id->calibr_fctr = hsp[i].cal_cst; id->pixel_xsize = fi->pixdim[1]; id->pixel_ysize = fi->pixdim[2]; bytes = id->width * id->height * MdcType2Bytes(id->type); if ( (id->buf=MdcGetImgBuffer(bytes)) == NULL ) { MdcFree(hsp); return("INW Bad malloc image buffer"); } /* read the image */ if (fread(id->buf,1,bytes,fi->ifp) != bytes) { err=MdcHandleTruncated(fi,i+1,MDC_YES); if (err != NULL) return(err); break; } } /* get the slice_width */ if (fi->number == 1) { /* unknown slice_width */ }else{ /* slice_width = relative translation from next slice */ fi->pixdim[3] = 0.0; for (i=0; inumber-1; i++) { if (hsp[i].trans > hsp[i+1].trans) { fi->image[i].slice_width = hsp[i].trans - hsp[i+1].trans; }else{ fi->image[i].slice_width = hsp[i+1].trans - hsp[i].trans; } fi->pixdim[3] += fi->image[i].slice_width; } /* last slice_width = previous slice_width */ fi->image[i].slice_width = fi->image[i-1].slice_width; fi->pixdim[3] += fi->image[i].slice_width; fi->pixdim[3] /= (float)fi->number; /* average */ } /* fill in orientation information */ fi->pat_slice_orient = MDC_SUPINE_HEADFIRST_TRANSAXIAL; /* default! */ strcpy(fi->pat_pos,MdcGetStrPatPos(fi->pat_slice_orient)); strcpy(fi->pat_orient,MdcGetStrPatOrient(fi->pat_slice_orient)); /* fill in the Acr/Nema variables */ for (i=0; inumber; i++) { id = &fi->image[i]; id->slice_spacing = id->slice_width; MdcFillImgPos(fi,i,i,0.0); MdcFillImgOrient(fi,i); } MdcFree(hsp); MdcCloseFile(fi->ifp); if (fi->truncated) return("INW Truncated image file"); return NULL; } int MdcWriteHeadStart(FILEINFO *fi) { MDC_INW_HEAD_START hs; memset(&hs,0,MDC_INW_HEAD_START_SIZE); hs.mark = MDC_INW_SIG; hs.version = MDC_INW_VERS_HIGH*256 + MDC_INW_VERS_LOW; hs.size_header = (Int16)( (MDC_INW_HEAD_SPEC_SIZE * fi->number) + MDC_INW_HEAD_START_SIZE + MDC_INW_HEAD_GEN_SIZE); hs.size_start = MDC_INW_HEAD_START_SIZE; hs.size_gen = MDC_INW_HEAD_GEN_SIZE; hs.size_spec = MDC_INW_HEAD_SPEC_SIZE; memcpy(hs.reserved,"MEDCON",6); MdcSWAP(hs.mark); MdcSWAP(hs.version); MdcSWAP(hs.size_header); MdcSWAP(hs.size_start); MdcSWAP(hs.size_gen); MdcSWAP(hs.size_spec); if (fwrite((char *)&hs,1,MDC_INW_HEAD_START_SIZE,fi->ofp) != MDC_INW_HEAD_START_SIZE) return(MDC_NO); return(MDC_YES); } int MdcWriteHeadGen(FILEINFO *fi) { MDC_INW_HEAD_GEN hg; memset(&hg,0,MDC_INW_HEAD_GEN_SIZE); hg.no = (Int16)fi->number; hg.sizeX = (Int16)fi->mwidth; hg.sizeY = (Int16)fi->mheight; hg.pixel_type = 2; /* only this supported */ hg.init_trans = 0; hg.dummy1 = 0; hg.time = 0; hg.decay_cst = 0.; hg.pixel_size = (fi->pixdim[1] + fi->pixdim[2]) / 2.; hg.max = (float) fi->qglmax; hg.min = (float) fi->qglmin; hg.scanner = 0; hg.reconstruction = '?'; hg.recon_version = 0; sprintf(hg.day,"%.11s ",MDC_DATE); sprintf(hg.reserved,"%.23s",MDC_PRGR); MdcMakeVAXfl(hg.decay_cst); MdcMakeVAXfl(hg.pixel_size); MdcMakeVAXfl(hg.max); MdcMakeVAXfl(hg.min); MdcSWAP(hg.no); MdcSWAP(hg.sizeX); MdcSWAP(hg.sizeY); MdcSWAP(hg.pixel_type); MdcSWAP(hg.init_trans); MdcSWAP(hg.dummy1); MdcSWAP(hg.time); MdcSWAP(hg.scanner); if (fwrite((char *)&hg,1,MDC_INW_HEAD_GEN_SIZE,fi->ofp) != MDC_INW_HEAD_GEN_SIZE) return(MDC_NO); return(MDC_YES); } int MdcSkipHeadSpecs(FILEINFO *fi) { Uint32 i; MDC_INW_HEAD_SPEC hsp; memset(&hsp,0,MDC_INW_HEAD_SPEC_SIZE); for (i=0; i < fi->number; i++) if (fwrite((char *)&hsp,1,MDC_INW_HEAD_SPEC_SIZE,fi->ofp) != MDC_INW_HEAD_SPEC_SIZE) return MDC_NO; return MDC_YES; } int MdcWriteHeadSpecs(FILEINFO *fi) { IMG_DATA *id; MDC_INW_HEAD_SPEC hsp; Uint32 img; fseek(fi->ofp,MDC_INW_HEAD_START_SIZE + MDC_INW_HEAD_GEN_SIZE,SEEK_SET); for (img=0; img < fi->number; img++) { memset(&hsp,0,MDC_INW_HEAD_SPEC_SIZE); id = &fi->image[img]; hsp.time = 0; if (id->rescaled) { hsp.max = (Int16)id->rescaled_max; hsp.min = (Int16)id->rescaled_min; hsp.cal_cst = id->rescaled_fctr; }else{ hsp.max = (Int16)id->max; hsp.min = (Int16)id->min; hsp.cal_cst = id->rescale_slope; } hsp.trans = img * (Int16)id->slice_width; #ifdef MDC_USE_SLICE_SPACING if (fi->number > 1) hsp.trans = img * (Int16)id->slice_spacing; #endif memcpy(hsp.reserved,MDC_INSTITUTION,6); MdcMakeVAXfl(hsp.cal_cst); MdcSWAP(hsp.time); MdcSWAP(hsp.trans); MdcSWAP(hsp.max); MdcSWAP(hsp.min); if (fwrite((char *)&hsp,1,MDC_INW_HEAD_SPEC_SIZE,fi->ofp) != MDC_INW_HEAD_SPEC_SIZE) return(MDC_NO); } return(MDC_YES); } char *MdcWriteINW(FILEINFO *fi) { IMG_DATA *id; double value; Uint32 i, p, size; Uint8 *buf, *maxbuf; int FREE=MDC_NO, type = BIT16_S; /* only supported type in version 1.0 */ MDC_FILE_ENDIAN = MDC_LITTLE_ENDIAN; /* always */ if (MDC_FORCE_INT != MDC_NO) { if (MDC_FORCE_INT != BIT16_S) { MdcPrntWarn("INW Only Int16 pixels supported"); } } if (XMDC_GUI == MDC_NO) { MdcDefaultName(fi,MDC_FRMT_INW,fi->ofname,fi->ifname); } if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_BEGIN,0.,"Writing INW:"); if (MDC_VERBOSE) MdcPrntMesg("INW Writing <%s> ...",fi->ofname); /* check for colored files */ if (fi->map == MDC_MAP_PRESENT) return("INW Colored files unsupported"); if (MDC_FILE_STDOUT == MDC_YES) { fi->ofp = stdout; }else{ if (MdcKeepFile(fi->ofname)) return("INW File exists!!"); if ( (fi->ofp=fopen(fi->ofname,"wb")) == NULL) return("INW Couldn't open file"); } if ( !MdcWriteHeadStart(fi) )return("INW Bad write HeadStart struct"); if ( !MdcWriteHeadGen(fi) ) return("INW Bad write HeadGen struct"); if ( !MdcSkipHeadSpecs(fi) ) return("INW Bad skipping HeadSpecs structs"); /* write the images */ for (i=0; inumber; i++) { if (MDC_PROGRESS) MdcProgress(MDC_PROGRESS_INCR,1./(float)fi->number,NULL); id = &fi->image[i]; if ((id->type != BIT16_S) || MDC_QUANTIFY || MDC_CALIBRATE) { buf = MdcGetImgBIT16_S(fi, i); FREE=MDC_YES; type=BIT16_S; }else{ buf = id->buf; FREE=MDC_NO; type=id->type; } if (buf == NULL) return("INW Bad malloc image buffer"); /* write images with uniform sizes */ if (fi->diff_size) { size = fi->mwidth * fi->mheight * MdcType2Bytes(type); maxbuf = MdcGetResizedImage(fi, buf, type, i); if (maxbuf == NULL) return("INW Bad malloc maxbuf"); if (FREE) MdcFree(buf); FREE = MDC_YES; }else{ size = id->width * id->height * MdcType2Bytes(type); maxbuf = buf; } for (p=0; p < size; p += MdcType2Bytes(type)) { value=MdcGetDoublePixel((Uint8 *)&maxbuf[p],type); MdcWriteDoublePixel(value,type,fi->ofp); } if (FREE) MdcFree(maxbuf); if (ferror(fi->ofp)) return("INW Bad images MdcFlush"); } if ( !MdcWriteHeadSpecs(fi) ) return("INW Bad write HeadSpecs structs"); MdcCheckQuantitation(fi); MdcCloseFile(fi->ofp); return NULL; }