tpcimage.h File Reference

Header file for libtpcimage. More...

#include "tpcclibConfig.h"
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "tpcextensions.h"
#include "tpcift.h"
#include "tpcisotope.h"
#include "tpctac.h"
#include "tpcecat.h"
#include "tpcdcm.h"
#include "tpcnifti.h"

Go to the source code of this file.

Data Structures
struct	IMG

Enumerations
enum	imgformat {   IMG_FORMAT_UNKNOWN , IMG_FORMAT_DICOM , IMG_FORMAT_E63 , IMG_FORMAT_E7 ,   IMG_FORMAT_E7_2D , IMG_FORMAT_POLARMAP , IMG_FORMAT_ANA , IMG_FORMAT_ANA_L ,   IMG_FORMAT_INTERFILE , IMG_FORMAT_NIFTI_1D , IMG_FORMAT_NIFTI_1S , IMG_FORMAT_NIFTI_2D ,   IMG_FORMAT_NIFTI_2S , IMG_FORMAT_MICROPET , IMG_FORMAT_FLAT , IMG_FORMAT_LAST }
enum	imgmodality {   IMG_MODALITY_UNKNOWN , IMG_MODALITY_PET , IMG_MODALITY_SPECT , IMG_MODALITY_CT ,   IMG_MODALITY_MRI , IMG_MODALITY_LAST }
enum	imgcontent {   IMG_CONTENT_UNKNOWN , IMG_CONTENT_IMAGE , IMG_CONTENT_RAW , IMG_CONTENT_ATTN ,   IMG_CONTENT_POLARMAP , IMG_CONTENT_LAST }

Functions
void	imgInit (IMG *image)
void	imgFree (IMG *image)
int	imgAllocate (IMG *img, const unsigned int dimz, const unsigned int dimy, const unsigned int dimx, const unsigned int dimt, TPCSTATUS *status)
char *	imgContentDescr (imgcontent c)
char *	imgModalityDescr (imgmodality c)
int	imgHasData (IMG *img)
int	imgHasTimes (IMG *img)
int	imgHasCounts (IMG *img)
int	imgHasWeights (IMG *img)
void	imgContents (IMG *img, FILE *fp)
unsigned long long	imgNaNs (IMG *img, int fix)
int	imgMinMax (IMG *img, float *minvalue, float *maxvalue)
int	imgXRange (IMG *img, double *xmin, double *xmax)
int	imgXUnitConvert (IMG *img, const int u)
int	imgCopyHeader (IMG *img1, IMG *img2)
int	imgFillOHeader (IMG *img, TPCSTATUS *status)
char *	imgFormatDescr (imgformat c)
char *	imgDefaultExtension (imgformat c)
int	imgRead (IMG *img, const char *fname, TPCSTATUS *status)
int	imgWrite (IMG *img, const char *fname, TPCSTATUS *status)
imgformat	imgFormatFromFName (const char *fname, TPCSTATUS *status)
imgformat	imgFormatIdentify (const char *s)
int	imgReadDICOM (IMG *img, const char *fname, TPCSTATUS *status)
int	imgWriteDICOM (IMG *img, const char *fname, TPCSTATUS *status)
int	niftiCreateFNames (const char *filename, char *hdrfile, char *imgfile, char *siffile, int fileformat)
int	imgReadNifti (IMG *img, const char *fname, TPCSTATUS *status)
int	imgWriteNifti (IMG *img, const char *fname, TPCSTATUS *status)
int	imgSetNiftiHeader (IMG *img, NIFTI_DSR *dsr, int verbose)
int	imgGetNiftiHeader (IMG *img, NIFTI_DSR *dsr, int verbose)
int	imgToSIF (IMG *img, TAC *sif, int copy_header, int copy_frames, int copy_counts, int verbose)
int	imgFromSIF (IMG *img, TAC *sif, int copy_header, int copy_frames, int copy_counts, int verbose)
int	imgCompareMatrixSize (IMG *d1, IMG *d2)
int	imgCompareUnit (IMG *d1, IMG *d2, TPCSTATUS *status)
int	imgCompareConc (IMG *d1, IMG *d2, const float test_abs, const float test_rel, TPCSTATUS *status)
int	imgCompareTimes (IMG *d1, IMG *d2, const float test_abs, const float test_rel, TPCSTATUS *status)

Detailed Description

Header file for libtpcimage.

Header file for image data processing library libtpcimage.

Author

Vesa Oikonen

Definition in file tpcimage.h.

Enumeration Type Documentation

imgcontent

enum imgcontent

Image data content codes for IMG data structure.

See also

img_content, image.c

Enumerator
IMG_CONTENT_UNKNOWN	Unknown data content.
IMG_CONTENT_IMAGE	Image data.
IMG_CONTENT_RAW	Sinogram, normalization.
IMG_CONTENT_ATTN	Attenuation data.
IMG_CONTENT_POLARMAP	Polarmap.
IMG_CONTENT_LAST	End of list.

Definition at line 70 of file tpcimage.h.

             {
  IMG_CONTENT_UNKNOWN,   
  IMG_CONTENT_IMAGE,     
  IMG_CONTENT_RAW,       
  IMG_CONTENT_ATTN,      
  IMG_CONTENT_POLARMAP,  
  IMG_CONTENT_LAST       
} imgcontent;

imgformat

enum imgformat

Image file format codes for IMG data structure.

See also

img_format, image.c

Enumerator
IMG_FORMAT_UNKNOWN	Unknown format.
IMG_FORMAT_DICOM	DICOM.
IMG_FORMAT_E63	ECAT 6.3.
IMG_FORMAT_E7	ECAT 7.* 3D.
IMG_FORMAT_E7_2D	ECAT 7.* 2D.
IMG_FORMAT_POLARMAP	ECAT polar map.
IMG_FORMAT_ANA	Analyze big-endian.
IMG_FORMAT_ANA_L	Analyze little-endian.
IMG_FORMAT_INTERFILE	Interfile.
IMG_FORMAT_NIFTI_1D	NIfTI-1 dual-file format.
IMG_FORMAT_NIFTI_1S	NIfTI-1 single-file format.
IMG_FORMAT_NIFTI_2D	NIfTI-2 dual-file format.
IMG_FORMAT_NIFTI_2S	NIfTI-2 single-file format.
IMG_FORMAT_MICROPET	MicroPET.
IMG_FORMAT_FLAT	Flat format.
IMG_FORMAT_LAST	End of list.

Definition at line 32 of file tpcimage.h.

             {
  IMG_FORMAT_UNKNOWN,     
  IMG_FORMAT_DICOM,       
  IMG_FORMAT_E63,         
  IMG_FORMAT_E7,          
  IMG_FORMAT_E7_2D,       
  IMG_FORMAT_POLARMAP,    
  IMG_FORMAT_ANA,         
  IMG_FORMAT_ANA_L,       
  IMG_FORMAT_INTERFILE,   
  IMG_FORMAT_NIFTI_1D,    
  IMG_FORMAT_NIFTI_1S,    
  IMG_FORMAT_NIFTI_2D,    
  IMG_FORMAT_NIFTI_2S,    
  IMG_FORMAT_MICROPET,    
  IMG_FORMAT_FLAT,        
  IMG_FORMAT_LAST         
} imgformat;

imgmodality

enum imgmodality

Image modality codes for IMG data structure.

See also

img_modality, image.c

Enumerator
IMG_MODALITY_UNKNOWN	Unknown modality.
IMG_MODALITY_PET	PET.
IMG_MODALITY_SPECT	SPECT.
IMG_MODALITY_CT	CT.
IMG_MODALITY_MRI	MRI.
IMG_MODALITY_LAST	End of list.

Definition at line 56 of file tpcimage.h.

             {
  IMG_MODALITY_UNKNOWN, 
  IMG_MODALITY_PET,     
  IMG_MODALITY_SPECT,   
  IMG_MODALITY_CT,      
  IMG_MODALITY_MRI,     
  IMG_MODALITY_LAST     
} imgmodality;

Function Documentation

imgAllocate()

int imgAllocate ( IMG * img, const unsigned int dimz, const unsigned int dimy, const unsigned int dimx, const unsigned int dimt, TPCSTATUS * status )

extern

Allocates memory for image data. Old contents are not saved.

See also

imgInit, imgFree, imgWrite, imgHasData, imgCopyHeader

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to initialized image structure; old contents are deleted.

Precondition

Initialize structure using imgInit() before this.

Parameters

dimz	Nr of image planes to allocate.
dimy	Nr of image rows to allocate.
dimx	Nr of image columns to allocate.
dimt	Nr of image time frames (samples) to allocate.
status	Pointer to status data; enter NULL if not needed.

Definition at line 126 of file image.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {
    printf("%s(img, %u, %u, %u, %u)\n", __func__, dimz, dimy, dimx, dimt); fflush(stdout);}
 
  if(img==NULL || dimz<1 || dimy<1 || dimx<1 || dimt<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
 
  /* Remove any previous contents */
  imgFree(img);
 
  /* Memory for frame times */
  int ret=0; {img->x1=calloc(dimt, sizeof(float)); if(img->x1==NULL) ret++;}
  if(!ret) {img->x2=calloc(dimt, sizeof(float)); if(img->x2==NULL) ret++;}
  if(!ret) {img->x=calloc(dimt, sizeof(float)); if(img->x==NULL) ret++;}
  if(ret) {
    imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
 
  /* Memory for weights */
  ret=0;   {img->weight=calloc(dimt, sizeof(float)); if(img->weight==NULL) ret++;}
  if(!ret) {img->prompts=calloc(dimt, sizeof(float)); if(img->prompts==NULL) ret++;}
  if(!ret) {img->randoms=calloc(dimt, sizeof(float)); if(img->randoms==NULL) ret++;}
  if(ret) {
    imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
 
  /* Memory for pixel matrix */
  ret=0;   {img->_t=(float*)malloc((size_t)dimz*dimy*dimx*dimt*sizeof(float)); if(img->_t==NULL) ret++;}
  if(!ret) {img->_x=(float**)malloc((size_t)dimz*dimy*dimx*sizeof(float*)); if(img->_x==NULL) ret++;}
  if(!ret) {img->_y=(float***)malloc((size_t)dimz*dimy*sizeof(float**)); if(img->_y==NULL) ret++;}
  if(!ret) {img->_z=(float****)malloc((size_t)dimz*sizeof(float***)); if(img->_z==NULL) ret++;}
  if(ret) {
    imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
  /* Set matrix data pointers */
  {
    float ***yptr, **xptr, *tptr;
    yptr=img->_y; xptr=img->_x; tptr=img->_t;
    for(unsigned int zi=0; zi<dimz; zi++) {
      img->_z[zi]=yptr;
      for(unsigned int yi=0; yi<dimy; yi++) {
        *yptr++=xptr;
        for(unsigned int xi=0; xi<dimx; xi++) {
          *xptr++=tptr; tptr+=dimt;
        }
      }
    }
  }
  img->m=img->_z;
  img->p=img->_t;
  /* Set all pixel values to zero */
  for(unsigned int zi=0; zi<dimz; zi++)
    for(unsigned int yi=0; yi<dimy; yi++)
      for(unsigned int xi=0; xi<dimx; xi++)
        for(unsigned int ti=0; ti<dimt; ti++)
          img->m[zi][yi][xi][ti]=(float)0.0;
  /* Set matrix dimensions */
  img->dimz=dimz; img->dimy=dimy; img->dimx=dimx; img->dimt=dimt;
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by imgReadDICOM(), and imgReadNifti().

imgCompareConc()

int imgCompareConc ( IMG * d1, IMG * d2, const float test_abs, const float test_rel, TPCSTATUS * status )

extern

Check whether IMG pixel concentrations are the same. Note that units are ignored here.

If either absolute or relative difference is below the limit, the test is reported as passed.

See also

tacCompareConc, imgCompareUnit, imgCompareTimes.

Returns

0 in case of match, and >0 if no match or error.

Author

Vesa Oikonen

Parameters

d1	Pointer to IMG structure.
d2	Pointer to IMG structure.
test_abs	Limit for accepted absolute difference; obligatory.
test_rel	Optional limit for accepted relative difference |2*(x1-x2)/(x1+x2)| ; set to negative value to not test this; in case of zero mean, this test is assumed to fail, but test for absolute difference may still pass.
status	Pointer to status data; enter NULL if not needed.

Definition at line 79 of file imagecomp.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>1) printf("%s()\n", __func__);
 
  /* Check that required data exists */
  if(d1==NULL || d2==NULL || d1->dimz<1 || d2->dimz<1 || d1->dimy<1 || d2->dimy<1
     || d1->dimx<1 || d2->dimx<1 || d1->dimt<1 || d2->dimt<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return 1;
  }
 
  /* Dimensions must match */
  if(d1->dimz!=d2->dimz || d1->dimy!=d2->dimy || d1->dimx!=d2->dimx || d1->dimt!=d2->dimt) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    if(verbose>0) printf("different IMG dimension.\n");
    return(2);
  }
 
  /* Compare */
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  for(int zi=0; zi<d1->dimz; zi++)
    for(int yi=0; yi<d1->dimy; yi++)
      for(int xi=0; xi<d1->dimx; xi++)
        for(int ti=0; ti<d1->dimt; ti++) {
          if(floatMatch(d1->m[zi][yi][xi][ti], d2->m[zi][yi][xi][ti], test_abs)==1)
            continue;
          if(test_rel>0.0 && floatMatchRel(d1->m[zi][yi][xi][ti], d2->m[zi][yi][xi][ti], test_rel)==1)
            continue;
          if(verbose>0) {
            float s;
            s=fabsf(d1->m[zi][yi][xi][ti] - d2->m[zi][yi][xi][ti]);
            printf("img1.m[%d][%d][%d][%d] := %g\n", zi, yi, xi, ti, d1->m[zi][yi][xi][ti]);
            printf("img2.m[%d][%d][%d][%d] := %g\n", zi, yi, xi, ti, d2->m[zi][yi][xi][ti]);
            printf("|diff| := %g\n", s);
            printf("diff_limit := %g\n", test_abs);
            if(test_rel>0.0) printf("rel_diff_limit := %g\n", test_rel);
          }
          return(10);
        }
  return(0);
}

imgCompareMatrixSize()

int imgCompareMatrixSize ( IMG * d1, IMG * d2 )

extern

Check whether two IMG data have the same matrix (x,y,z) size.

Returns

0 in case of match, 1 otherwise.

See also

imgCompareConc, imgCompareTimes

Author

Vesa Oikonen

Parameters

d1	Pointer to IMG structure.
d2	Pointer to IMG structure.

Definition at line 18 of file imagecomp.c.

  {
  if(d1==NULL || d2==NULL) return(1);
  if(d1->dimz!=d2->dimz) return(1);
  if(d1->dimy!=d2->dimy) return(1);
  if(d1->dimx!=d2->dimx) return(1);
  return(0);
}

imgCompareTimes()

int imgCompareTimes ( IMG * d1, IMG * d2, const float test_abs, const float test_rel, TPCSTATUS * status )

extern

Check whether IMG frame times (x values) are the same. Note that units are ignored here.

If either absolute or relative difference is below the limit, the test is reported as passed.

See also

tacCompareTimes, imgCompareConc, imgCompareUnit, tacimgXMatch.

Returns

0 in case of match, and >0 if no match or error.

Author

Vesa Oikonen

Parameters

d1	Pointer to IMG structure.
d2	Pointer to IMG structure.
test_abs	Limit for accepted absolute difference.
test_rel	Optional limit for accepted relative difference |2*(x1-x2)/(x1+x2)| ; set to negative value to not test this; in case of zero mean, this test is assumed to fail, but test for absolute difference may still pass.
status	Pointer to status data; enter NULL if not needed.

Definition at line 144 of file imagecomp.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>1) printf("%s()\n", __func__);
 
  /* Check that required data exists */
  if(d1==NULL || d2==NULL || d1->dimt<1 || d2->dimt<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return 1;
  }
 
  /* Frame nr must match */
  if(d1->dimt!=d2->dimt) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    if(verbose>0) printf("different frame nr.\n");
    return(3);
  }
  
  /* Compare */
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  for(int fi=0; fi<d1->dimt; fi++) {
    if(floatMatch(d1->x[fi], d2->x[fi], test_abs)==1)
      continue;
    if(floatMatch(d1->x1[fi], d2->x1[fi], test_abs)==1)
      continue;
    if(floatMatch(d1->x2[fi], d2->x2[fi], test_abs)==1)
      continue;
    if(test_rel>0.0 && 
       floatMatchRel(d1->x[fi], d2->x[fi], test_rel)==1 &&
       floatMatchRel(d1->x1[fi], d2->x1[fi], test_rel)==1 &&
       floatMatchRel(d1->x2[fi], d2->x2[fi], test_rel)==1)
      continue;
    if(verbose>0) {
      printf("img1.x1[%d] := %g\n", fi, d1->x1[fi]);
      printf("img2.x1[%d] := %g\n", fi, d2->x1[fi]);
      printf("|diff| := %g\n", fabsf(d1->x1[fi]-d2->x1[fi]));
      printf("img1.x2[%d] := %g\n", fi, d1->x2[fi]);
      printf("img2.x2[%d] := %g\n", fi, d2->x2[fi]);
      printf("|diff| := %g\n", fabsf(d1->x2[fi]-d2->x2[fi]));
      printf("img1.x[%d] := %g\n", fi, d1->x[fi]);
      printf("img2.x[%d] := %g\n", fi, d2->x[fi]);
      printf("|diff| := %g\n", fabsf(d1->x[fi]-d2->x[fi]));
      printf("diff_limit := %g\n", test_abs);
      if(test_rel>0.0) printf("rel_diff_limit := %g\n", test_rel);
    }
    return(11);
  }
  return(0);
}

imgCompareUnit()

int imgCompareUnit ( IMG * d1, IMG * d2, TPCSTATUS * status )

extern

Check whether time and concentration units are the same in two IMG data.

Returns

0 in case of match, 1 if y (concentration) unit is not matching, 2 if x (time) unit is not matching, and 3 if neither is matching.

See also

tacCompareUnit, imgCompareConc, imgCompareTimes

Author

Vesa Oikonen

Parameters

d1	Pointer to IMG structure.
d2	Pointer to IMG structure.
status	Pointer to status data; enter NULL if not needed

Definition at line 39 of file imagecomp.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>1) printf("%s()\n", __func__);
 
  /* Check that required data exists */
  if(d1==NULL || d2==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(3);
  }
 
  int ret=0;
  if(d1->cunit!=d2->cunit) ret+=1;
  if(d1->tunit!=d2->tunit) ret+=2;
  if(ret>0 && verbose>0) {
    printf("img1.cunit := %s\n", unitName(d1->cunit));
    printf("img2.cunit := %s\n", unitName(d2->cunit));
    printf("img1.tunit := %s\n", unitName(d1->tunit));
    printf("img2.tunit := %s\n", unitName(d2->tunit));
  }
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(ret);
}

imgContentDescr()

char * imgContentDescr ( imgcontent c )

extern

Return pointer to content type of the image data.

Returns

pointer to the data content type string.

See also

imgdata

Parameters

c	IMG content code

Definition at line 51 of file image.c.

  {
  if(c<IMG_CONTENT_UNKNOWN || c>=IMG_CONTENT_LAST) return NULL;
  return (char*)img_content[c];
}

Referenced by imgContents(), and imgFillOHeader().

imgContents()

void imgContents ( IMG * img, FILE * fp )

extern

Print general information on the contents of the IMG data structure.

See also

imgInit, imgRead, imgHasData, imgHasTimes

Parameters

img	Pointer to IMG structure.
fp	File pointer for output; usually stdout.

Definition at line 293 of file image.c.

  {
  if(fp==NULL) return;
  if(img==NULL) {fprintf(fp, "IMG pointer is NULL\n"); fflush(fp); return;}
 
  fprintf(fp, "modality: %s\n", imgModalityDescr(img->modality));
  fprintf(fp, "content: %s\n", imgContentDescr(img->content));
  fprintf(fp, "format: %s\n", imgFormatDescr(img->format));
  if(img->oformat!=IMG_FORMAT_UNKNOWN)
    fprintf(fp, "output_format: %s\n", imgFormatDescr(img->oformat));
 
  fprintf(fp, "hasData := %d\n", imgHasData(img));
  fprintf(fp, "hasTimes := %d\n", imgHasTimes(img));
  fprintf(fp, "hasCounts := %d\n", imgHasCounts(img));
 
  if(img->scanStart[0]) fprintf(fp, "scanStart: %s\n", img->scanStart);
 
  fprintf(fp, "dimensions: %u x %u x %u x %u\n", img->dimx, img->dimy, img->dimz, img->dimt);
  if(img->sizex!=0.0 || img->sizey!=0.0 || img->sizez!=0.0)
    fprintf(fp, "voxel_size: %g x %g x %g\n", img->sizex, img->sizey, img->sizez);
  if(img->gapx!=0.0 || img->gapy!=0.0 || img->gapz!=0.0)
    fprintf(fp, "voxel_gaps: %g x %g x %g\n", img->gapx, img->gapy, img->gapz);
 
  fprintf(fp, "qform := %d\n", img->xform[0]);
  fprintf(fp, "sform := %d\n", img->xform[1]);
  if(floatNonzeroes(img->quatern, 6)>0) {
    fprintf(fp, "quatern_b := %g\n", img->quatern[0]);
    fprintf(fp, "quatern_c := %g\n", img->quatern[1]);
    fprintf(fp, "quatern_d := %g\n", img->quatern[2]);
    fprintf(fp, "quatern_x_shift := %g\n", img->quatern[3]);
    fprintf(fp, "quatern_y_shift := %g\n", img->quatern[4]);
    fprintf(fp, "quatern_z_shift := %g\n", img->quatern[5]);
  }
  if(floatNonzeroes(img->srow, 12)>0) {
    for(int i=0; i<4; i++) fprintf(fp, "srow_x[%d] := %g\n", i, img->srow[i]);
    for(int i=0; i<4; i++) fprintf(fp, "srow_y[%d] := %g\n", i, img->srow[4+i]);
    for(int i=0; i<4; i++) fprintf(fp, "srow_z[%d] := %g\n", i, img->srow[8+i]);
  }
  if(floatNonzeroes(img->iop, 6)>0) {
    fprintf(fp, "image_orientation_patient := %g", img->iop[0]);
    for(int i=1; i<6; i++) fprintf(fp, ", %g", img->iop[i]);
    fprintf(fp, "\n");
  }
  if(floatNonzeroes(img->ipp, 3)>0) {
    fprintf(fp, "image_position_patient := %g", img->ipp[0]);
    for(int i=1; i<3; i++) fprintf(fp, ", %g", img->ipp[i]);
    fprintf(fp, "\n");
  }
  if(floatNonzeroes(img->mt, 3)>0) {
    fprintf(fp, "image_matrix_transformation_parameters := %g", img->mt[0]);
    for(int i=1; i<12; i++) fprintf(fp, ", %g", img->mt[i]);
    fprintf(fp, "\n");
  }
 
  if(img->isot!=ISOTOPE_UNKNOWN) fprintf(fp, "isotope: %s\n", isotopeName(img->isot));
  if(img->decayCorrection==DECAY_CORRECTED) fprintf(fp, "physical_decay: corrected\n");
  else if(img->decayCorrection==DECAY_NOTCORRECTED) fprintf(fp, "physical_decay: not corrected\n");
 
  if(imgHasWeights(img)) fprintf(fp, "weighting: no\n");
  else fprintf(fp, "weighting: yes\n");
 
  if(img->cunit!=UNIT_UNKNOWN) fprintf(fp, "pixel_unit := %s\n", unitName(img->cunit));
  if(img->tunit!=UNIT_UNKNOWN) fprintf(fp, "time_unit := %s\n", unitName(img->tunit));
 
  fprintf(fp, "original header: %d items\n", img->ih.keyNr);
 
  fflush(fp);
  return;
}

imgCopyHeader()

int imgCopyHeader ( IMG * img1, IMG * img2 )

extern

Copy information fields from one IMG to another. Image dimensions, pixel data, and frame information are not copied.

See also

imgAllocate, imgInit, imgFree

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img1	Pointer to source image structure.
img2	Pointer to initiated target image structure.

Definition at line 18 of file imageheader.c.

  {
  if(img1==NULL || img2==NULL) return(TPCERROR_FAIL);
  if(!imgHasData(img1) || !imgHasData(img2)) return(TPCERROR_FAIL);
  
  strlcpy(img2->studyNr, img1->studyNr, MAX_STUDYNR_LEN+1);
  img2->isot=img1->isot;
  img2->decayCorrection=img1->decayCorrection;
  strlcpy(img2->scanStart, img1->scanStart, 20);
  img2->content=img1->content;
  img2->modality=img1->modality;
  img2->format=img1->format;
  img2->oformat=img1->oformat;
  img2->sizex=img1->sizex;
  img2->sizey=img1->sizey;
  img2->sizez=img1->sizez;
  img2->gapx=img1->gapx;
  img2->gapy=img1->gapy;
  img2->gapz=img1->gapz;
  img2->xform[0]=img1->xform[0];
  img2->xform[1]=img1->xform[1];
  for(int i=0; i<6; i++) img2->quatern[i]=img1->quatern[i];
  for(int i=0; i<12; i++) img2->srow[i]=img1->srow[i];
  for(int i=0; i<6; i++) img2->iop[i]=img1->iop[i];
  for(int i=0; i<3; i++) img2->ipp[i]=img1->ipp[i];
  for(int i=0; i<12; i++) img2->mt[i]=img1->mt[i];
  img2->cunit=img1->cunit;
  img2->tunit=img1->tunit;
 
  return(TPCERROR_OK);
}

imgDefaultExtension()

char * imgDefaultExtension ( imgformat c )

extern

Return pointer to default image file name extension based on the format code.

Remarks

No extension for most formats.

Returns

pointer to the file name extension starting with '.', or empty string.

See also

imgformat

Parameters

c	IMG format code

Definition at line 67 of file imageio.c.

  {
  if(c<IMG_FORMAT_UNKNOWN || c>=IMG_FORMAT_LAST) return NULL;
  return (char*)img_fn_ext[c];
}

imgFillOHeader()

int imgFillOHeader ( IMG * img, TPCSTATUS * status )

extern

Fill/edit IMG output header.

Remarks

Stub function.

See also

imgInit, imgWrite, imgFree

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to image structure.
status	Pointer to status data; enter NULL if not needed.

Definition at line 61 of file imageheader.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(img)\n", __func__); fflush(stdout);}
 
  if(img==NULL || !imgHasData(img)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(TPCERROR_NO_DATA);
  }
 
  /* If output header is empty then copy input header;
     if file formats match, copy directly, otherwise convert contents */
  if(img->oh.keyNr==0 && img->ih.keyNr>0) {
    if(img->format==img->oformat) {
      if(verbose>1) printf("copying header\n");
      if(iftDuplicate(&img->ih, &img->oh)!=TPCERROR_OK) {
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_HEADER);
        return(TPCERROR_INVALID_HEADER);
      }
    } else {
      if(verbose>1) printf("converting header\n");
 
    }
  }
 
  /* Set fields from IMG structure */
 
  iftDeleteKey(&img->oh, "dimx");
  if(iftPutUInt(&img->oh, "dimx", img->dimz, 0, NULL)!=0) return(TPCERROR_INVALID_HEADER);
  iftDeleteKey(&img->oh, "dimy");
  if(iftPutUInt(&img->oh, "dimy", img->dimz, 0, NULL)!=0) return(TPCERROR_INVALID_HEADER);
  iftDeleteKey(&img->oh, "dimz");
  if(iftPutUInt(&img->oh, "dimz", img->dimz, 0, NULL)!=0) return(TPCERROR_INVALID_HEADER);
  iftDeleteKey(&img->oh, "dimt");
  if(iftPutUInt(&img->oh, "dimt", img->dimz, 0, NULL)!=0) return(TPCERROR_INVALID_HEADER);
 
  if(img->modality!=IMG_MODALITY_UNKNOWN) {
    iftDeleteKey(&img->oh, "modality");
    if(iftPut(&img->oh, "modality", imgModalityDescr(img->modality), 0, NULL)!=0) 
      return(TPCERROR_INVALID_HEADER);
  }
  if(img->content!=IMG_CONTENT_UNKNOWN) {
    iftDeleteKey(&img->oh, "content");
    if(iftPut(&img->oh, "content", imgContentDescr(img->content), 0, NULL)!=0) 
      return(TPCERROR_INVALID_HEADER);
  }
  if(img->format!=IMG_FORMAT_UNKNOWN) {
    iftDeleteKey(&img->oh, "format");
    if(iftPut(&img->oh, "format", imgFormatDescr(img->format), 0, NULL)!=0) 
      return(TPCERROR_INVALID_HEADER);
  }
  if(img->isot!=ISOTOPE_UNKNOWN) {
    iftDeleteKey(&img->oh, "isotope");
    if(iftPut(&img->oh, "isotope", isotopeName(img->isot), 0, NULL)!=0)
      return(TPCERROR_INVALID_HEADER);
  }
  if(img->decayCorrection!=DECAY_UNKNOWN) {
    iftDeleteKey(&img->oh, "decay_correction");
    char buf[8];
    if(img->decayCorrection==DECAY_CORRECTED) strcpy(buf, "y"); else strcpy(buf, "n");
    if(iftPut(&img->oh, "decay_correction", buf, 0, NULL)!=0)
      return(TPCERROR_INVALID_HEADER);
  }
  if(img->cunit!=UNIT_UNKNOWN) {
    iftDeleteKey(&img->oh, "pixel_unit");
    if(iftPut(&img->oh, "pixel_unit", unitName(img->cunit), 0, NULL)!=0)
      return(TPCERROR_INVALID_HEADER);
  }
  if(img->tunit!=UNIT_UNKNOWN) {
    iftDeleteKey(&img->oh, "time_unit");
    if(iftPut(&img->oh, "time_unit", unitName(img->tunit), 0, NULL)!=0)
      return(TPCERROR_INVALID_HEADER);
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by imgWrite().

imgFormatDescr()

char * imgFormatDescr ( imgformat c )

extern

Return pointer to image file format description with the format code.

Returns

pointer to the image file format string.

See also

imgformat

Parameters

c	IMG format code

Definition at line 36 of file imageio.c.

  {
  if(c<IMG_FORMAT_UNKNOWN || c>=IMG_FORMAT_LAST) return NULL;
  return (char*)img_format[c];
}

Referenced by imgContents(), imgFillOHeader(), imgReadNifti(), and imgWrite().

imgFormatFromFName()

imgformat imgFormatFromFName ( const char * fname, TPCSTATUS * status )

extern

Determine IMG file format from file name.

If file name is actually name of an existing path, then DICOM is assumed. Non-existing path can be identified as path only if it ends with '/'.

Note that NIfTI dual file format, Analyze, and microPET files can not be separated by file naming.

Returns

enum imgformat (IMG_FORMAT_DICOM, ...), and IMG_FORMAT_UNKNOWN, when not certain.

See also

imgFormatIdentify, imgRead, imgWrite, micropetExists, niftiExists, anaExists

Parameters

fname	Name of file that is used to determine format.
status	Pointer to status data; enter NULL if not needed.

Definition at line 394 of file imageio.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(%s)\n", __func__, fname); fflush(stdout);}
 
  /* Empty name? */
  if(fname==NULL || strnlen(fname, 2)<1) return(IMG_FORMAT_UNKNOWN);
 
  /* Path as file name? */
  if(pathExist(fname)) {
    if(verbose>1) printf("%s is an existing path\n", fname);
    return(IMG_FORMAT_DICOM);
  }
  if(fname[strlen(fname)-1]=='/' || fname[strlen(fname)-1]=='\\') {
    if(verbose>1) printf("%s is a non-existing path\n", fname);
    return(IMG_FORMAT_DICOM);
  }
 
  /* Get the extension(s) */
  char *exts=filenameGetExtensions(fname);
  if(exts==NULL) {
    if(verbose>1) printf("file name has no extension\n");
    /* Probably NIfTI, Analyze, or MicroPET, but we cannot be sure */
    return(IMG_FORMAT_UNKNOWN);
  }
  if(verbose>1) printf("file name extension: '%s'\n", exts);
 
  /* Check certain common extensions */
  if(strcasecmp(exts, ".dcm")==0) return(IMG_FORMAT_DICOM);
  if(strcasecmp(exts, ".dicom")==0) return(IMG_FORMAT_DICOM);
 
  if(strcasecmp(exts, ".nii")==0) return(IMG_FORMAT_NIFTI_1S);
 
  if(strcasecmp(exts, ".v")==0) return(IMG_FORMAT_E7);
  if(strcasecmp(exts, ".s")==0) return(IMG_FORMAT_E7);
  if(strcasecmp(exts, ".n")==0) return(IMG_FORMAT_E7);
  if(strcasecmp(exts, ".a")==0) return(IMG_FORMAT_E7);
  if(strcasecmp(exts, ".scn")==0) return(IMG_FORMAT_E63);
  if(strcasecmp(exts, ".nrm")==0) return(IMG_FORMAT_E63);
  if(strcasecmp(exts, ".atn")==0) return(IMG_FORMAT_E63);
  if(strcasestr(exts, "polarmap")!=NULL) return(IMG_FORMAT_POLARMAP);
 
  if(strcasecmp(exts, ".if")==0) return(IMG_FORMAT_INTERFILE);
  if(strcasestr(exts, ".i.hdr")!=NULL) return(IMG_FORMAT_INTERFILE);
  if(strcasestr(exts, ".i.img")!=NULL) return(IMG_FORMAT_INTERFILE);
  if(strcasecmp(exts, ".i")==0) return(IMG_FORMAT_INTERFILE);
 
  if(strcasecmp(exts, ".img")==0) return(IMG_FORMAT_E63);
 
  if(strcasecmp(exts, ".pet.img")==0) return(IMG_FORMAT_MICROPET);
  if(strcasecmp(exts, ".pet.img.hdr")==0) return(IMG_FORMAT_MICROPET);
  if(strcasecmp(exts, ".ct.img")==0) return(IMG_FORMAT_MICROPET);
  if(strcasecmp(exts, ".ct.img.hdr")==0) return(IMG_FORMAT_MICROPET);
 
  /* Check certain in-house extensions */
  if(strcasecmp(exts, ".bin")==0) return(IMG_FORMAT_FLAT);
  if(strcasecmp(exts, ".mif")==0) return(IMG_FORMAT_FLAT);
 
  return(IMG_FORMAT_UNKNOWN);
}

Referenced by imgWrite().

imgFormatIdentify()

imgformat imgFormatIdentify ( const char * s )

extern

Identify the string representation of the IMG file format.

Returns

enum imgformat, or 0 (enum IMG_FORMAT_UNKNOWN) if not identified.

See also

imgFormatFromFName

Author

Vesa Oikonen

Parameters

s	IMG format as a string

Definition at line 466 of file imageio.c.

  {
  if(s==NULL || strlen(s)<1) return IMG_FORMAT_UNKNOWN;
  /* Try if string can be found directly in the table */
  for(int i=0; i<IMG_FORMAT_LAST; i++) {
    if(strcasecmp(img_format[i], s)==0) return i;
  }
  /* Format string is not following TPC standard, lets try something else */
  if(     strcasecmp(s, "bin")==0)              return IMG_FORMAT_FLAT;
  else if(strcasecmp(s, "v")==0)                return IMG_FORMAT_E7;
  else if(strcasecmp(s, "dcm")==0)              return IMG_FORMAT_DICOM;
  else if(strcasecmp(s, "ana")==0)              return IMG_FORMAT_ANA_L;
  else if(strcasecmp(s, "nii")==0)              return IMG_FORMAT_NIFTI_1S;
  else if(strcasecmp(s, "nifti")==0)            return IMG_FORMAT_NIFTI_1S;
 
  return IMG_FORMAT_UNKNOWN;
}

imgFree()

void imgFree ( IMG * img )

extern

Free memory allocated for IMG. All data is cleared.

See also

imgInit, imgAllocate, imgRead

Parameters

img	Pointer to IMG structure.

Definition at line 107 of file image.c.

  {
  if(img==NULL) return;
  iftFree(&img->ih);
  iftFree(&img->oh);
  free(img->weight); free(img->prompts); free(img->randoms);
  free(img->x1); free(img->x2); free(img->x);
  free(img->_t); free(img->_z); free(img->_x); free(img->_y);
  imgInit(img);
}

Referenced by imgAllocate(), imgRead(), imgReadDICOM(), and imgReadNifti().

imgFromSIF()

int imgFromSIF ( IMG * img, TAC * sif, int copy_header, int copy_frames, int copy_counts, int verbose )

extern

Copy information from SIF into IMG.

Todo

Add tests.

See also

imgToSIF, imgInit, tacInit, imgWrite, imgRead, tacWriteSIF

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to IMG struct into which content is copied from SIF.
sif	Pointer to SIF, which actually is stored as a TAC format.
copy_header	Select whether header contents are copied (1) or not copied (0) from SIF.
copy_frames	Select whether frame times are copied (1) or not copied (0) from SIF.
copy_counts	Select whether counts are copied (1) or not copied (0) from SIF.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout

Definition at line 126 of file imagesif.c.

  {
  if(verbose>0) {
    printf("%s(img, sif, %d, %d, %d, ...)\n", __func__, copy_header, copy_frames, copy_counts);
    fflush(stdout);
  }
 
  if(sif==NULL || img==NULL) return(TPCERROR_FAIL);
 
  if(copy_header) {
    if(verbose>1) printf("  copying header.\n");
    int ret=tacGetHeaderScanstarttime(&sif->h, img->scanStart, NULL);
    if(ret!=TPCERROR_OK && verbose>0) fprintf(stderr, "Warning: cannot read scan start time in SIF\n");
    img->isot=tacGetIsotope(sif);
    if(img->isot==ISOTOPE_UNKNOWN && verbose>0) fprintf(stderr, "Warning: cannot read isotope in SIF\n");
    ret=tacGetHeaderStudynr(&sif->h, img->studyNr, NULL);
  }
 
  if(copy_frames) {
    if(verbose>1) printf("  copying frame times.\n");
    if(!sif->isframe || sif->sampleNr!=img->dimt) return(TPCERROR_INVALID_FORMAT);
    for(int fi=0; fi<img->dimt; fi++) {
      img->x1[fi]=sif->x1[fi];
      img->x2[fi]=sif->x2[fi];
      img->x[fi]=0.5*(img->x1[fi]+img->x2[fi]);
    }
    img->tunit=sif->tunit;
  }
 
  if(copy_counts) {
    if(verbose>1) printf("  copying count data.\n");
    if(sif->format!=TAC_FORMAT_SIF || !sif->isframe || sif->sampleNr!=img->dimt)
      return(TPCERROR_INVALID_FORMAT);
    if(sif->tacNr<4) return(TPCERROR_INVALID_FORMAT);
    for(int fi=0; fi<img->dimt; fi++) {
      img->prompts[fi]=sif->c[0].y[fi];
      img->randoms[fi]=sif->c[1].y[fi];
    }
    img->weighting=WEIGHTING_OFF;
    if(sif->weighting!=WEIGHTING_OFF && sif->weighting!=WEIGHTING_UNKNOWN) {
      for(int fi=0; fi<img->dimt; fi++)
        img->weight[fi]=sif->w[fi];
      img->weighting=sif->weighting;
    }
  }
 
  return(TPCERROR_OK);
}

Referenced by imgReadNifti().

imgGetNiftiHeader()

int imgGetNiftiHeader ( IMG * img, NIFTI_DSR * dsr, int verbose )

extern

Copy header information from NIfTI into IMG.

Todo

Add tests.

See also

imgSetNiftiHeader, niftiExists, imgReadNifti

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to IMG structure to be filled.
dsr	Pointer to NIfTI header structure from which information is read.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 794 of file imagenii.c.

  {
  if(verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
  
  /* Check the input */
  if(dsr==NULL || img==NULL) return(TPCERROR_FAIL);
  if(dsr->n!=1 && dsr->n!=2) return(TPCERROR_FAIL);
 
  /* NIfTI file format */
  if(dsr->n==1) {
    if(strcmp(dsr->h1.magic, "ni1")==0) img->format=IMG_FORMAT_NIFTI_1D;
    else if(strcmp(dsr->h1.magic, "n+1")==0) img->format=IMG_FORMAT_NIFTI_1S;
    else return(TPCERROR_INVALID_FORMAT);
  } else if (dsr->n==2) {
    if(strcmp(dsr->h2.magic, "ni2")==0) img->format=IMG_FORMAT_NIFTI_2D;
    else if(strcmp(dsr->h2.magic, "n+2")==0) img->format=IMG_FORMAT_NIFTI_2S;
    else return(TPCERROR_INVALID_FORMAT);
  } else {
    return(TPCERROR_FAIL);
  }
 
  /* Get the image dimensions */
  int dimNr;
  if(dsr->n==1) dimNr=dsr->h1.dim[0]; else dimNr=dsr->h2.dim[0];
  if(dimNr<2 || dimNr>8) return(TPCERROR_INVALID_HEADER);
  if(dimNr>4) {
    if(verbose>0) fprintf(stderr, "Error: NIfTI image dimension %d is not supported\n", dimNr);
    return(TPCERROR_UNSUPPORTED);
  }
 
  long long dimx, dimy, dimz=1, dimt=1;
  if(dsr->n==1) {
    dimx=dsr->h1.dim[1];
    dimy=dsr->h1.dim[2];
    if(dimNr>2) {dimz=dsr->h1.dim[3]; if(dimNr>3) dimt=dsr->h1.dim[4];}
  } else {
    dimx=dsr->h2.dim[1];
    dimy=dsr->h2.dim[2];
    if(dimNr>2) {dimz=dsr->h2.dim[3]; if(dimNr>3) dimt=dsr->h2.dim[4];}
  }
  long long pxlNr=dimx*dimy*dimz;
  if(pxlNr<1 || dimt<0) {
    if(verbose>0) fprintf(stderr, "Error: invalid NIfTI image dimensions.\n");
    return(TPCERROR_INVALID_HEADER);
  }
  if(dimt==0) dimt=1;
  img->dimx=dimx; img->dimy=dimy; img->dimz=dimz; img->dimt=dimt;
 
  /* Pixel x,y,z sizes, converting units if necessary */
  int xyzt_units=0;
  if(dsr->n==1) xyzt_units=dsr->h1.xyzt_units; else xyzt_units=dsr->h2.xyzt_units;
  float f=1.0;
  if(xyzt_units & NIFTI_UNITS_METER) f=1000.;
  else if(xyzt_units & NIFTI_UNITS_MICRON) f=0.001;
  else if(xyzt_units & NIFTI_UNITS_MM) f=1.0;
  if(verbose>2) printf("pixel size conversion factor := %g\n", f);
  if(dsr->n==1) {
    img->sizex=dsr->h1.pixdim[1]; img->sizey=dsr->h1.pixdim[2]; img->sizez=dsr->h1.pixdim[3];
  } else {
    img->sizex=dsr->h2.pixdim[1]; img->sizey=dsr->h2.pixdim[2]; img->sizez=dsr->h2.pixdim[3];
  }
  img->sizex*=f; img->sizey*=f; img->sizez*=f;
 
  /* Orientation, quaternion, and transformation parameters */
  if(dsr->n==1) {
    img->xform[0]=dsr->h1.qform_code;
    img->xform[1]=dsr->h1.sform_code;
    img->quatern[0]=dsr->h1.quatern_b;
    img->quatern[1]=dsr->h1.quatern_c;
    img->quatern[2]=dsr->h1.quatern_d;
    img->quatern[3]=dsr->h1.qoffset_x;
    img->quatern[4]=dsr->h1.qoffset_y;
    img->quatern[5]=dsr->h1.qoffset_z;
    for(int i=0; i<4; i++) img->srow[i]=dsr->h1.srow_x[i];
    for(int i=0; i<4; i++) img->srow[4+i]=dsr->h1.srow_y[i];
    for(int i=0; i<4; i++) img->srow[8+i]=dsr->h1.srow_z[i];
  } else {
    img->xform[0]=dsr->h2.qform_code;
    img->xform[1]=dsr->h2.sform_code;
    img->quatern[0]=dsr->h2.quatern_b;
    img->quatern[1]=dsr->h2.quatern_c;
    img->quatern[2]=dsr->h2.quatern_d;
    img->quatern[3]=dsr->h2.qoffset_x;
    img->quatern[4]=dsr->h2.qoffset_y;
    img->quatern[5]=dsr->h2.qoffset_z;
    for(int i=0; i<4; i++) img->srow[i]=dsr->h2.srow_x[i];
    for(int i=0; i<4; i++) img->srow[4+i]=dsr->h2.srow_y[i];
    for(int i=0; i<4; i++) img->srow[8+i]=dsr->h2.srow_z[i];
  }
 
  /* Assumptions etc */
  img->content=IMG_CONTENT_IMAGE;
  strcpy(img->studyNr, "");
  img->decayCorrection=DECAY_CORRECTED;
  img->cunit=UNIT_UNKNOWN;
 
  /* Check if unit can be found in description field */
  char *cptr;
  if(dsr->n==1) cptr=dsr->h1.descrip; else cptr=dsr->h2.descrip;
  img->cunit=unitIdentify(cptr);
 
  return(TPCERROR_OK);
}

Referenced by imgReadNifti().

imgHasCounts()

int imgHasCounts ( IMG * img )

extern

Verify that IMG contains prompts and randoms.

See also

imgHasData, imgHasTimes, imgHasWeights, imgInit, imgFree, imgRead, imgWrite, imgContents

Returns

0, if no prompts or randoms, 1 if prompts exist, 2 if randoms exist, and 3 if both prompts and randoms exist.

Parameters

img	Pointer to image structure.

Definition at line 251 of file image.c.

  {
  if(img==NULL || img->dimt<1) return(0);
  int p=0, r=0;
  /* If image has only one frame, then any value > 0 is fine */
  if(img->dimt==1) {
    if(img->prompts[0]>0.000000001) p=1;
    if(img->randoms[0]>0.000000001) r=2;
    return(p+r);
  }
  /* When more frames, we check that frames have different count level */
  for(unsigned short int i=1; i<img->dimt; i++) {
    if(fabsf(img->prompts[i]-img->prompts[i-1])>0.001) p=1;
    if(fabsf(img->randoms[i]-img->randoms[i-1])>0.001) r=2;
    if((p+r)>2) break;
  }
  return(p+r);
}

Referenced by imgContents(), and imgToSIF().

imgHasData()

int imgHasData ( IMG * img )

extern

Verify that image data actually contains pixel data.

See also

imgHasTimes, imgNaNs, imgRead, imgWrite, imgContents

Returns

0, if no data, 1 if contains data.

Parameters

img	Pointer to image structure.

Definition at line 218 of file image.c.

  {
  if(img==NULL) return(0);
  if(img->dimz<1 || img->dimy<1 || img->dimx<1 || img->dimz<1) return(0);
  if(img->m==NULL || img->p==NULL) return(0);
  if(img->x1==NULL || img->x2==NULL) return(0);
  return(1);
}

Referenced by imgContents(), imgCopyHeader(), imgFillOHeader(), imgHasTimes(), imgMinMax(), imgSetNiftiHeader(), imgWrite(), imgWriteNifti(), and imgXRange().

imgHasTimes()

int imgHasTimes ( IMG * img )

extern

Verify that IMG contains frame times.

See also

imgHasData, imgInit, imgFree, imgRead, imgWrite, imgContents

Returns

0, if no frame times, 1 if frame times exist.

Parameters

img	Pointer to image structure.

Definition at line 235 of file image.c.

  {
  if(imgHasData(img)==0) return(0);
  for(int i=0; i<img->dimt; i++) if(img->x2[i]>0.00000001) return(1);
  return(0);
}

Referenced by imgContents(), imgToSIF(), and imgWriteNifti().

imgHasWeights()

int imgHasWeights ( IMG * img )

extern

Check if image contains weights as indicated by the 'weighting' field.

See also

imgHasCounts, imgHasData, imgHasTimes, imgInit, tacIsWeighted

Returns

Returns 1 if weighted, and 0 if not or in case of an error.

Parameters

img	Pointer to image structure.

Definition at line 278 of file image.c.

  {
  if(img==NULL || img->dimt<1) return(0);
  if(img->weighting==WEIGHTING_UNKNOWN || img->weighting==WEIGHTING_OFF) return(0);
  if(img->weighting<WEIGHTING_LAST) return(1);
  return(0);
}

Referenced by imgContents().

imgInit()

void imgInit ( IMG * img )

extern

Initiate the IMG structure before any use.

See also

imgFree, imgAllocate, imgRead

Parameters

img	Pointer to IMG structure.

Definition at line 64 of file image.c.

  {
  if(img==NULL) return;
  img->dimt=img->dimx=img->dimy=img->dimz=0;
  iftInit(&img->ih);
  iftInit(&img->oh);
  img->content=IMG_CONTENT_UNKNOWN;
  img->modality=IMG_MODALITY_UNKNOWN;
  img->format=img->oformat=IMG_FORMAT_UNKNOWN;
 
  for(int i=0; i<20; i++) img->scanStart[i]=(char)0;
  for(int i=0; i<MAX_STUDYNR_LEN+1; i++) img->studyNr[i]=(char)0;
 
 
  img->sizex=img->sizey=img->sizez=0.0;
  img->gapx=img->gapy=img->gapz=0.0;
  img->xform[0]=img->xform[1]=0;
  for(int i=0; i<6; i++) img->quatern[i]=0.0;
  for(int i=0; i<12; i++) img->srow[i]=0.0;
  for(int i=0; i<6; i++) img->iop[i]=0.0;
  for(int i=0; i<3; i++) img->ipp[i]=0.0;
  for(int i=0; i<12; i++) img->mt[i]=0.0;
 
  img->isot=ISOTOPE_UNKNOWN;
  img->decayCorrection=DECAY_UNKNOWN;
 
  img->weighting=WEIGHTING_UNKNOWN;
  img->weight=img->prompts=img->randoms=NULL;
 
  img->cunit=img->tunit=UNIT_UNKNOWN;
 
  img->x1=img->x2=img->x=NULL;
  img->m=(float****)0; img->p=0;
  img->_z=(float****)0; img->_y=(float***)0; img->_x=(float**)0; img->_t=(float*)0;
}

Referenced by imgFree().

imgMinMax()

int imgMinMax ( IMG * img, float * minvalue, float * maxvalue )

extern

Searches the min and max pixel value in the image data.

See also

imgXRange, imgInit, imgFree, imgHasData

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to IMG structure from where min and max pixel values are searched.
minvalue	Pointer to min pixel value; Enter NULL if not needed.
maxvalue	Pointer to max pixel value; Enter NULL if not needed.

Definition at line 396 of file image.c.

  {
  if(img==NULL || !imgHasData(img)) return(TPCERROR_FAIL);
  float f1, f2;
  f1=f2=nanf("");
  if(minvalue!=NULL) *minvalue=f1;
  if(maxvalue!=NULL) *maxvalue=f2;
  unsigned long long i, n=img->dimt*img->dimz*img->dimy*img->dimx;
  // search first valid number
  for(i=0; i<n; i++) if(isfinite(img->p[i])) break;
  if(i==n) return(TPCERROR_MISSING_DATA);
  f1=f2=img->p[i++];
  for(; i<n; i++) {
    if(img->p[i]>f2) f2=img->p[i];
    else if(img->p[i]<f1) f1=img->p[i];
  }
  if(minvalue!=NULL) *minvalue=f1;
  if(maxvalue!=NULL) *maxvalue=f2;
  return(TPCERROR_OK);
}

Referenced by imgSetNiftiHeader().

imgModalityDescr()

char * imgModalityDescr ( imgmodality c )

extern

Return pointer to modality of the image data.

Returns

pointer to the modality string.

See also

imgmodality

Parameters

c	IMG modality code

Definition at line 27 of file image.c.

  {
  if(c<IMG_MODALITY_UNKNOWN || c>=IMG_MODALITY_LAST) return NULL;
  return (char*)img_modality[c];
}

Referenced by imgContents(), and imgFillOHeader().

imgNaNs()

unsigned long long imgNaNs ( IMG * img, int fix )

extern

Searches the image data for missing pixel values, optionally setting those to zero.

See also

imgHasData

Returns

the number of missing pixel values.

Parameters

img	Pointer to IMG structure.
fix	Set (1) or do not set (0) missing pixels to zero.

Definition at line 373 of file image.c.

  {
  if(img==NULL) return(0);
  unsigned long long i, n=0, pxlNr=img->dimt*img->dimz*img->dimy*img->dimx;
  for(i=0; i<pxlNr; i++) {
    if(!isfinite(img->p[i])) {
      n++;
      if(fix!=0) img->p[i]=0.0;
    }
  }
  return(n);
}

imgRead()

int imgRead ( IMG * img, const char * fname, TPCSTATUS * status )

extern

Read image file into IMG structure.

Remarks

Stub function, mostly.

See also

imgInit, imgWrite, imgFree

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to initialized image structure; old contents are deleted.

Precondition

Initialize structure using imgInit() before this.

Parameters

fname	Pointer to the file name, possibly path in case of DICOM.
status	Pointer to status data; enter NULL if not needed.

Definition at line 82 of file imageio.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(img, %s)\n", __func__, fname); fflush(stdout);}
 
  if(img==NULL || strnlen(fname, 2)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
 
  /* Delete any old contents in data structure */
  imgFree(img);
 
  /* What kind of data are we supposed to read? */
  if(pathExist(fname)) {
 
    if(verbose>2) printf("%s is a path\n", fname);
    unsigned short int fileNr=pathFileNr(fname);
    if(verbose>2) printf("%s contains %u file(s)\n", fname, fileNr);
    if(fileNr==0) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_FILE);
      return(TPCERROR_NO_FILE);
    }
    /* Check if this is a DICOM folder */
    IFT fl; iftInit(&fl);
    int ret=dcmFileList(fname, &fl, status);
    if(ret==TPCERROR_OK) {
      if(verbose>2) printf("%s is DICOM folder.\n", fname);
      img->format=IMG_FORMAT_DICOM;
    } else {
      if(verbose>2) printf("%s is not a DICOM folder.\n", fname);
      fileNr=pathFileList(fname, &fl); // just read the file list
    }
    if(fileNr==0) {
      iftFree(&fl);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_FILE);
      return(TPCERROR_NO_FILE);
    }
    /* It would be too risky to just open random file from the folder */
    iftFree(&fl);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return(TPCERROR_UNSUPPORTED);
 
  } else if(fileExist(fname)) {
 
    if(verbose>2) printf("%s is a file\n", fname);
    /* Check if this is DICOM file */
    if(dcmVerifyMagic(fname, NULL)) {
      /*if(verbose>2)*/ printf("%s is in DICOM format.\n", fname);
      img->format=IMG_FORMAT_DICOM;
      /* This may be a single-file DICOM, or user just gave the name of one file in DICOM folder */
      IFT fl; iftInit(&fl);
      int ret=dcmFileList(fname, &fl, status);
      if(ret!=TPCERROR_OK) {iftFree(&fl); return(ret);}
      if(fl.keyNr==1) {
        if(verbose>2) printf("  %s is a single DICOM file\n", fname);
        iftFree(&fl);
        ret=imgReadDICOM(img, fname, status);
        if(verbose>2) printf("  imgReadDICOM() := %s\n", errorMsg(status->error));
        return(ret);
      } else {
        if(verbose>2) printf("  %s is one of %d DICOM files\n", fname, fl.keyNr);
 
        iftFree(&fl);
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
        return(TPCERROR_UNSUPPORTED);
      }
    }
    /* Check if we have NIfTI file, which may be in single file format,
       or dual file format which has similar names as microPET and Analyze */
    if(niftiExists(fname, NULL, NULL, NULL, NULL, status)==1) {
      if(verbose>2) printf("%s is in NIfTI format.\n", fname);
      /* Read NIfTI image */
      return(imgReadNifti(img, fname, status));
      /*statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
      return(TPCERROR_UNSUPPORTED);*/
    }
    /* Check if we have Analyze file, which may be in dual file format which 
       has similar names as microPET and dual-file NIfTI */
    if(anaExists(fname, NULL, NULL, NULL, NULL, status)==1) {
      /* Read Analyze image */
      /*if(verbose>2)*/ printf("%s is in Analyze format.\n", fname);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
      return(TPCERROR_UNSUPPORTED);
    }
 
  } else { // Not a folder, neither an existing file
 
    /* Check if this is base name of an image, without extension(s) */
    /* Check if we have NIfTI file, which may be in single file format,
       or dual file format which has similar names as microPET and Analyze */
    if(niftiExists(fname, NULL, NULL, NULL, NULL, status)==1) {
      /* Read NIfTI image */
      /*if(verbose>2)*/ printf("%s is in NIfTI format.\n", fname);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
      return(TPCERROR_UNSUPPORTED);
    }
    /* Check if we have Analyze file, which may be in dual file format which 
       has similar names as microPET and dual-file NIfTI */
    if(anaExists(fname, NULL, NULL, NULL, NULL, status)==1) {
      /* Read Analyze image */
      /*if(verbose>2)*/ printf("%s is in Analyze format.\n", fname);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
      return(TPCERROR_UNSUPPORTED);
    }
    /* Check for microPET */
 
 
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_FILE);
    return(TPCERROR_NO_FILE);
  }
 
 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

imgReadDICOM()

int imgReadDICOM ( IMG * img, const char * fname, TPCSTATUS * status )

extern

Read DICOM PET image into IMG data structure.

Remarks

Only few header fields are read.

Todo

Add tests.

See also

imgInit, imgRead, imgWrite, imgFree, imgWriteDICOM

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to image structure. Any previous contents removed.
fname	Pointer to the DICOM file name. Path not accepted.
status	Pointer to status data; enter NULL if not needed.

Definition at line 23 of file imagedcm.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(..., %s, ...)\n", __func__, fname); fflush(stdout);}
 
  if(img==NULL || strnlen(fname, 2)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
 
  /* Delete any old contents in data structure */
  imgFree(img);
 
  /* Read DICOM file */
  DCMFILE dcm; dcmfileInit(&dcm);
  int ret=dcmFileRead(fname, &dcm, 0, status);
  if(ret!=TPCERROR_OK) {dcmfileFree(&dcm); return(ret);}
 
//verbose=20;
 
  /* Get modality */
  imgmodality modality=IMG_MODALITY_UNKNOWN;
  {
    DCMTAG tag; dcmTagSet(&tag, 0x0008, 0x0060); DCMITEM *iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL) {
      if(verbose>0) {fprintf(stderr, "Error: cannot find modality.\n"); fflush(stderr);}
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    char *buf=dcmValueString(iptr); 
    if(verbose>12) printf("Modality := '%s'\n", buf);
    if(strncmp(buf, "PT", 2)==0) modality=IMG_MODALITY_PET;
    else if(strncmp(buf, "CT", 2)==0) modality=IMG_MODALITY_CT;
    else if(strncmp(buf, "MR", 2)==0) modality=IMG_MODALITY_MRI;
    else if(strncmp(buf, "NM", 2)==0) modality=IMG_MODALITY_SPECT;
    else if(verbose>0) fprintf(stderr, "Warning: modality '%s' not supported.\n", buf);
    free(buf);
  }
 
  /* Get the number of frames and planes */
  if(verbose>10) {printf("reading matrix dimensions\n"); fflush(stdout);}
  unsigned short int imgdim[4];
  ret=dcmImgDim(&dcm, imgdim, verbose-2);
  if(ret!=0 || imgdim[0]<=0 || imgdim[1]<=0 || imgdim[2]<=0 || imgdim[3]<=0) {
    if(verbose>0) {fprintf(stderr, "Error: cannot find matrix dimensions.\n"); fflush(stderr);}
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
    dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
  }
  unsigned short int dimx, dimy, sliceNr, frameNr;
  dimx=imgdim[0]; dimy=imgdim[1]; sliceNr=imgdim[2]; frameNr=imgdim[3];
  if(verbose>11) 
    printf("dimx := %u\ndimy := %u\ndimz := %u\ndimt := %u\n", dimx, dimy, sliceNr, frameNr); 
 
  /* 
   *  Allocate memory for IMG data
   */
  ret=imgAllocate(img, sliceNr, dimy, dimx, frameNr, status);
  if(ret!=TPCERROR_OK) {dcmfileFree(&dcm); return(ret);}
 
 
  /* Get the status of decay correction and isotope */
  ret=dcmImgIsotope(&dcm, &img->isot, &img->decayCorrection, verbose-2);
  if(ret!=0 && verbose>0) {
    fprintf(stderr, "Warning: decay correction and isotope unknown.\n"); fflush(stderr);
  }
  if(verbose>11) {
    printf("decayCorrection := %s\n", decayDescr(img->decayCorrection));
    printf("isotope := %s\n", isotopeName(img->isot));
    fflush(stdout);
  }
 
  /* If decay is corrected, read time of decay correction; otherwise read start time */
  char zeroDateTime[32]; zeroDateTime[0]=(char)0;
  if(img->decayCorrection>=DECAY_CORRECTED) {
    DCMTAG tag; dcmTagSet(&tag, 0x0018, 0x9701);
    DCMITEM *iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) dcmDT2intl(iptr->rd, zeroDateTime);
  } else {
    /* Get Acquisition DateTime */
    DCMTAG tag; dcmTagSet(&tag, 0x0008, 0x002A);
    DCMITEM *iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) dcmDT2intl(iptr->rd, zeroDateTime);
  }
  if(!zeroDateTime[0]) {
    if(verbose>0) fprintf(stderr, "Error: missing Acquisition Date and Time.\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
    imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
  }
  if(verbose>14) printf("zeroDateTime := %s\n", zeroDateTime);
 
  /* Get frame times */
  if(modality==IMG_MODALITY_PET || modality==IMG_MODALITY_SPECT) {
    /* Initiate frame times to NaN */
    for(unsigned short int i=0; i<img->dimt; i++) img->x1[i]=img->x2[i]=img->x[i]=nanf("");
 
    /* Find Per Frame Functional Groups Sequence */
    DCMTAG tag; dcmTagSet(&tag, 0x5200, 0x9230);
    DCMITEM *iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL) {
      if(verbose>0) fprintf(stderr, "Error: Per Frame Functional Groups Sequence not found.\n");
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    /* Loop through all Frame Content Sequences under it */
    dcmTagSet(&tag, 0x0020, 0x9111);
    DCMTAG tagstart; dcmTagSet(&tagstart, 0x0018, 0x9074);
    DCMTAG tagdur; dcmTagSet(&tagdur, 0x0018, 0x9220);
    DCMTAG tagfr; dcmTagSet(&tagfr, 0x0020, 0x9128);
    DCMITEM *fptr=iptr->child_item;
    while(fptr!=NULL) {
      /* Find Frame Content Sequence */
      iptr=dcmFindDownTag(fptr, 0, &tag, 0); if(iptr==NULL) break;
      if(verbose>20) printf(" found Frame Content Sequence\n");
      DCMITEM *jptr;
      /* Find Temporal Position Index */
      jptr=dcmFindDownTag(iptr->child_item, 0, &tagfr, 0);
      if(jptr==NULL) {fptr=fptr->next_item; continue;}
      if(verbose>30) dcmitemPrint(jptr);
      unsigned short int frameIndex=(unsigned short int)(dcmitemGetInt(jptr)-1);
      /* Find Frame Acquisition DateTime */
      jptr=dcmFindDownTag(iptr->child_item, 0, &tagstart, 0);
      char facqDateTime[32]; facqDateTime[0]=(char)0;
      if(jptr!=NULL) dcmDT2intl(jptr->rd, facqDateTime);
      if(!facqDateTime[0]) {
        if(verbose>0) fprintf(stderr, "Error: missing Frame Acquisition DateTime.\n");
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
        imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
      }
      if(verbose>20) printf("facqDateTime := %s\n", facqDateTime); 
      /* Calculate the Frame start time */
      float t1=(float)strDateTimeDifference(facqDateTime, zeroDateTime);
      /* Set frame start time */
      if(frameIndex<img->dimt) {
        if(isnan(img->x1[frameIndex])) {
          /* Save frame start time, since there is no previous value for this frame */
          img->x1[frameIndex]=t1;
          if(verbose>29) printf("t1[%u]=%g\n", frameIndex, t1);
        } else {
          /* Check whether previous frame start time is the same as current */
          if(!floatMatch(img->x1[frameIndex], t1, 0.001)) {
            // not the same, this probably is whole-body study with more than one bed pos
            if(verbose>1) printf("  different t1=%g\n", t1);
          }
        }
      }
 
      /* Find Frame Acquisition Duration */
      jptr=dcmFindDownTag(iptr->child_item, 0, &tagdur, 0);
      if(jptr==NULL) {fptr=fptr->next_item; continue;}
      if(verbose>20) dcmitemPrint(jptr);
      if(frameIndex<img->dimt) {
        img->x2[frameIndex]=img->x1[frameIndex]+(float)0.001*dcmitemGetReal(jptr);
        img->x[frameIndex]=0.5*(img->x1[frameIndex]+img->x2[frameIndex]);
      }
      /* prepare for the next frame */
      fptr=iptr->next_item;
    }
  }
 
  /* Get pixel sizes and image position */
  {
    double pxlsize[3];
    ret=dcmImgPxlsize(&dcm, pxlsize, verbose-2);
    if(ret!=0) {
      if(verbose>0) fprintf(stderr, "Error: missing Pixel size.\n");
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    img->sizex=(float)pxlsize[0];
    img->sizey=(float)pxlsize[1];
    img->sizez=(float)pxlsize[2];
    if(verbose>13) printf("  pixel_size := %g x %g x %g\n", img->sizex, img->sizey, img->sizez);
 
    double ipp[3];
    int ret1=dcmImgPos(&dcm, ipp, verbose-2);
    if(ret1==0) {
      for(int i=0; i<3; i++) img->ipp[i]=ipp[i];
    } else {
      if(verbose>0) fprintf(stderr, "Warning: missing Image Position (Patient).\n");
    }
    double iop[6];
    int ret2=dcmImgOrient(&dcm, iop, verbose-2);
    if(ret2==0) {
      for(int i=0; i<6; i++) img->iop[i]=iop[i];
    } else {
      if(verbose>0) fprintf(stderr, "Warning: missing Image Orientation (Patient).\n");
    }
    if(ret1==0 && ret2==0) { // calculate xform into IMG
      double xform[16];
      ret=dcmImgXform(iop, pxlsize, ipp, xform, verbose-2);
      if(ret==0) {
        img->xform[0]=1;
        for(int i=0; i<12; i++) img->srow[i]=(float)xform[i];
        double quatern[3], qoffset[3];
        if(dcmXformToQuatern(xform, quatern, qoffset, verbose-2)==0) {
          for(int i=0; i<3; i++) img->quatern[i]=(float)quatern[i];
          for(int i=0; i<3; i++) img->quatern[3+i]=(float)qoffset[i];
        }
      }
    }
  }
 
 
  /* Get Accession Number and write it in IMG as studyNr */
  {
    img->studyNr[0]=(char)0;
    DCMITEM *iptr; DCMTAG tag; dcmTagSet(&tag, 0x0008, 0x0050);
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) {
      char *buf=dcmValueString(iptr); if(verbose>12) printf("Accession Number := '%s'\n", buf);
      if(buf!=NULL && strlen(buf)>1) strncpyClean(img->studyNr, buf, MAX_STUDYNR_LEN+1);
      free(buf);
    }
    if(!strcasecmp(img->studyNr, "empty")) img->studyNr[0]=(char)0;
  }
  /* If not successful, then try Patient ID */
  if(!img->studyNr[0]) {
    DCMITEM *iptr; DCMTAG tag; dcmTagSet(&tag, 0x0010, 0x0020);
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr!=NULL) {
      char *buf=dcmValueString(iptr); if(verbose>12) printf("Patient ID := '%s'\n", buf);
      if(buf!=NULL && strlen(buf)>1) strncpyClean(img->studyNr, buf, MAX_STUDYNR_LEN+1);
      free(buf);
    }
  } 
 
 
  /* Get pixel data scaling factors */
  double rescale_intercept=0.0, rescale_slope=0.0;
  {
    DCMTAG tag; DCMITEM *iptr, *jptr;
    dcmTagSet(&tag, 0x0040, 0x9096);
    iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL) {
      if(verbose>0) fprintf(stderr, "Error: missing Real World Value Mapping Sequence.\n");
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_KEY);
      imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_KEY);
    }
    dcmTagSet(&tag, 0x0040, 0x9224);
    jptr=dcmFindDownTag(iptr->child_item, 0, &tag, 0);
    if(jptr!=NULL) rescale_intercept=dcmitemGetReal(jptr);
    dcmTagSet(&tag, 0x0040, 0x9225);
    jptr=dcmFindDownTag(iptr->child_item, 0, &tag, 0);
    if(jptr!=NULL) rescale_slope=dcmitemGetReal(jptr);
    if(rescale_slope<=0.0 && verbose>0) {
      fprintf(stderr, "Warning: rescale_slope := %g\n", rescale_slope); fflush(stderr);
    } else if(verbose>14) {
      printf("rescale_slope := %g\n", rescale_slope);
      printf("rescale_intercept := %g\n", rescale_intercept); // always 0 for PET
    }
    /* Unit */
    dcmTagSet(&tag, 0x0040, 0x9210);
    jptr=dcmFindDownTag(iptr->child_item, 0, &tag, 0);
    if(jptr!=NULL) {
      char *buf=dcmValueString(jptr); if(verbose>14) printf("  buf := '%s'\n", buf);
      strClean(buf); img->cunit=unitIdentify(buf);
      free(buf);
    }
    if(verbose>14) printf("image_unit := %s\n", unitName(img->cunit));
  }
 
  /* Fill IMG pixel values */
  {
    DCMTAG tag; dcmTagSet(&tag, 0x7fe0, 0x0010);
    DCMITEM *iptr=dcmFindTag(dcm.item, 0, &tag, 0);
    if(iptr==NULL || iptr->rd==NULL) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
      imgFree(img); dcmfileFree(&dcm); return(TPCERROR_NO_DATA);
    }
    unsigned long long bi=0, byteNr=iptr->vl;
    if(verbose>13) printf("  pixel data byteNr := %llu\n", byteNr);
 
    char *cptr=iptr->rd;
    short int s;
    for(unsigned int ti=0; ti<img->dimt; ti++) {
      for(unsigned int zi=0; zi<img->dimz; zi++)
        for(unsigned int yi=0; yi<img->dimy; yi++)
          for(unsigned int xi=0; xi<img->dimx; xi++) {
            if(bi>=byteNr) {printf("missing pixel data\n"); fflush(stdout); break;}
            memcpy(&s, cptr, 2); cptr+=2; bi+=2;
            img->m[zi][yi][xi][ti]=(float)(rescale_slope*(double)s + rescale_intercept);
          }
    }
  }
 
 
  dcmfileFree(&dcm);
 
  /* Set IMG header fields */
  img->modality=modality;
  img->format=IMG_FORMAT_DICOM;
  img->content=IMG_CONTENT_IMAGE;
  strcpy(img->scanStart, zeroDateTime);
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by imgRead().

imgReadNifti()

int imgReadNifti ( IMG * img, const char * fname, TPCSTATUS * status )

extern

Read NIfTI into IMG data structure.

Note

NIfTI may have missing pixel values (NaNs), remove those with imgNaNs(IMG, 1) when necessary.

Todo

Add tests.

See also

imgInit, imgRead, imgWrite, imgFree, imgNaNs, niftiCreateFNames

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to image structure. Any previous contents removed.
fname	Pointer to the database name. Path alone not accepted.
status	Pointer to status data; enter NULL if not needed.

Definition at line 72 of file imagenii.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(..., %s, ...)\n", __func__, fname); fflush(stdout);}
 
  if(img==NULL || strnlen(fname, 2)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
  /* Delete any old contents in data structure */
  imgFree(img);
 
  /* Get NIfTI filenames, check existence, and read header */
  if(verbose>1) printf("checking that %s exists and is NIfTI\n", fname);
  char hdrfile[FILENAME_MAX], imgfile[FILENAME_MAX], siffile[FILENAME_MAX];
  NIFTI_DSR dsr;
  if(!niftiExists(fname, hdrfile, imgfile, siffile, &dsr, status)) {
    if(status!=NULL) return(status->error);
    else return(TPCERROR_FAIL);
  }
  if(verbose>1) printf("NIfTI header read from %s\n", hdrfile);
 
  /* Get data type and check that it is currently supported */
  if(verbose>2) printf("  verifying datatype\n");
  int datatype=0;
  if(dsr.n==1) datatype=dsr.h1.datatype; else datatype=dsr.h2.datatype;
  switch(datatype) {
    case NIFTI_DT_UNSIGNED_CHAR:
    case NIFTI_DT_SIGNED_SHORT:
    case NIFTI_DT_SIGNED_INT:
    case NIFTI_DT_FLOAT:
    case NIFTI_DT_DOUBLE:
    case NIFTI_DT_SIGNED_CHAR:
    case NIFTI_DT_UNSIGNED_SHORT:
    case NIFTI_DT_UNSIGNED_INT:
    case NIFTI_DT_LONG_LONG:
    case NIFTI_DT_UNSIGNED_LONG_LONG:
      break;
    case NIFTI_DT_UNKNOWN:
    case NIFTI_DT_BINARY:
    case NIFTI_DT_COMPLEX:
    case NIFTI_DT_RGB:
    case NIFTI_DT_ALL:
    case NIFTI_DT_LONG_DOUBLE:
    case NIFTI_DT_DOUBLE_PAIR:
    case NIFTI_DT_LONG_DOUBLE_PAIR:
    case NIFTI_DT_RGBA:
    default:
      if(verbose>0) printf("Error: currently unsupported NIfTI datatype %d\n", datatype);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
      return(TPCERROR_UNSUPPORTED);
  }
  int bitpix=0;
  if(dsr.n==1) bitpix=dsr.h1.bitpix; else bitpix=dsr.h2.bitpix;
  if(verbose>2) printf("  bits per pixel := %d\n", bitpix);
  if(bitpix==0) { // Carimas Nifti Writer does not set bitpix
    if(datatype==NIFTI_DT_UNSIGNED_SHORT) bitpix=16;
  }
  if(bitpix<8 || (bitpix%8)!=0) { // We don't support bit data
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return(TPCERROR_UNSUPPORTED);
  }
  if(bitpix>64) { // We don't currently support 128 bit data or larger
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return(TPCERROR_UNSUPPORTED);
  }
  int bytepix=bitpix/8; if(verbose>2) printf("  bytes per pixel := %d\n", bytepix);
 
 
  /* Allocate memory for 4D image */
  {
  int ret=imgGetNiftiHeader(img, &dsr, verbose-1); // Get the image dimensions
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  if(verbose>1) printf("allocating memory for 4D image\n");
  ret=imgAllocate(img, img->dimz, img->dimy, img->dimx, img->dimt, status);
  if(ret) return(ret);
  }
 
  /* Copy header contents into IMG structure */
  {
  int ret=imgGetNiftiHeader(img, &dsr, verbose-1);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  if(verbose>2) printf("format: %s\n", imgFormatDescr(img->format));
  }
 
  /* Open file */
  if(verbose>1) printf("reading pixel data in %s\n", imgfile);
  FILE *fp=fopen(imgfile, "rb");
  if(fp==NULL) {
    imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_OPEN);
    return(TPCERROR_CANNOT_OPEN);
  }
 
 
  /* Set initial read position.
     Get the image data start location from header, in case of single file format */
  long long start_pos=0;
  {
    long long int s=0;
    if(dsr.n==1) s=(int)dsr.h1.vox_offset; else s=(int)dsr.h2.vox_offset;
    if(s<0) start_pos=-s; else start_pos=s;
    /* Check that start location is not earlier than it can by definition; it can be later, though */
    if(img->format==IMG_FORMAT_NIFTI_1S && s<NIFTI1_HEADER_SIZE) s=NIFTI1_HEADER_SIZE;
    else if(img->format==IMG_FORMAT_NIFTI_2S && s<NIFTI2_HEADER_SIZE) s=NIFTI2_HEADER_SIZE;
    start_pos=s;
  }
  if(verbose>2) printf("  image_start_pos := %llu\n", start_pos);
  /* Seek the data position, jumping over possible header */
  if(start_pos>0) {
    fseeko(fp, start_pos, SEEK_SET);
    if(ftello(fp)!=start_pos) {
      fclose(fp); imgFree(img);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FORMAT);
      return(TPCERROR_INVALID_FORMAT);
    }
  }
 
  /* Allocate memory for one frame of the binary data */
  long long pxlNr=img->dimx*img->dimy*img->dimz;
  if(verbose>2) printf("  pixels/frame := %llu\n", pxlNr);
  long long rawSize=pxlNr*(bitpix/8);
  if(verbose>1) printf("  raw bytes per frame := %lld\n", rawSize);
  if(verbose>1) printf("  allocating memory for raw binary data\n");
  unsigned char *buf=(unsigned char*)malloc(rawSize);
  if(buf==NULL) {
    fclose(fp); imgFree(img);
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    return(TPCERROR_OUT_OF_MEMORY);
  }
 
  /* Read the data, frame-by-frame */
  if(verbose>1) printf("  reading binary data\n");
  int little=endianLittle(); // Are we on little endian platform?
  int byteconv=0;
  if(little!=dsr.byte_order) {
    byteconv=1;
    if(verbose>1) printf("  byte conversion needed\n");
  }
  /* Get scaling factors */
  float scl_slope, scl_inter;
  if(dsr.n==1) {scl_slope=dsr.h1.scl_slope; scl_inter=dsr.h1.scl_inter;}
  else {scl_slope=dsr.h2.scl_slope; scl_inter=dsr.h2.scl_inter;}
  if(scl_slope==0.0) scl_slope=1.0;
 
  for(int ti=0; ti<img->dimt; ti++) {
    if(verbose>3) printf("  frame %d\n", 1+ti);
    if(fread(buf, rawSize, 1, fp) < 1) {
      fclose(fp); imgFree(img); free(buf);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_READ);
      return(TPCERROR_CANNOT_READ);
    }
    /* Convert byte order if necessary */
    if(byteconv) switch(bitpix) {
      case 8: /* no conversion needed */ break;
      case 16: swap16ip(buf, pxlNr); break;
      case 32: swap32ip(buf, pxlNr); break;
      case 64: swap64ip(buf, pxlNr); break;
      default: /* others not supported */ break;
    }
    /* Copy data to float pixel values */
    int ret=0;
    unsigned char *bptr=buf;
    switch(datatype) {
      case NIFTI_DT_UNSIGNED_CHAR:
      case NIFTI_DT_SIGNED_CHAR:
        if(bitpix!=8) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)(*bptr);
              bptr++;
            }
        break;
      case NIFTI_DT_UNSIGNED_SHORT:
        if(bitpix!=16) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              uint16_t a;
              memcpy(&a, bptr, 2);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_SIGNED_SHORT:
        if(bitpix!=16) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              int16_t a;
              memcpy(&a, bptr, 2);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_SIGNED_INT:
        if(bitpix!=32) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              int32_t a;
              memcpy(&a, bptr, 4);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_UNSIGNED_INT:
        if(bitpix!=32) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              uint32_t a;
              memcpy(&a, bptr, 4);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_LONG_LONG:
        if(bitpix!=64) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              int64_t a;
              memcpy(&a, bptr, 8);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_UNSIGNED_LONG_LONG:
        if(bitpix!=64) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              uint64_t a;
              memcpy(&a, bptr, 8);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_FLOAT:
        if(bitpix!=32) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              float a;
              memcpy(&a, bptr, 4);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*a;
              bptr+=bytepix;
            }
        break;
      case NIFTI_DT_DOUBLE:
        if(bitpix!=64) {ret=1; break;}
        for(int zi=0; zi<img->dimz; zi++)
          for(int yi=0; yi<img->dimy; yi++)
            for(int xi=0; xi<img->dimx; xi++) {
              double a;
              memcpy(&a, bptr, 8);
              img->m[zi][yi][xi][ti]=scl_inter+scl_slope*(float)a;
              bptr+=bytepix;
            }
        break;
      default:
        ret=2;
    }
    if(ret) {
      fclose(fp); imgFree(img); free(buf);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
      return(TPCERROR_UNSUPPORTED);
    }
  } // next frame
  fclose(fp); free(buf);
 
 
  /* Read SIF, if available */
  if(siffile[0]) {
    if(verbose>0) {printf("reading SIF %s\n", siffile); fflush(stdout);}
    TAC sif; tacInit(&sif);
    int ret=tacRead(&sif, siffile, status);
    if(ret!=TPCERROR_OK) {imgFree(img); return(ret);}
    ret=imgFromSIF(img, &sif, 1, 1, 0, verbose-1);
    tacFree(&sif);
    if(ret!=TPCERROR_OK) {
      imgFree(img);
      statusSet(status, __func__, __FILE__, __LINE__, ret);
      return(ret);
    }   
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by imgRead().

imgSetNiftiHeader()

int imgSetNiftiHeader ( IMG * img, NIFTI_DSR * dsr, int verbose )

extern

Copy header information in IMG structure into NIfTI header structure.

Header contents will depend on the output format as specified in IMG->oformat, which can be set to IMG_FORMAT_NIFTI_1D, IMG_FORMAT_NIFTI_1S, IMG_FORMAT_NIFTI_2D, or IMG_FORMAT_NIFTI_2S. NIFTI_DSR->n will tell whether NIfTI-1 or NIfTI-2 header was filled.

Todo

Add tests.

See also

imgGetNiftiHeader, imgWriteNifti, niftiWriteHeader

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to IMG structure from which header information is read.
dsr	Pointer to NIfTI header structure to be filled; any previous contents are removed.
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 576 of file imagenii.c.

  {
  if(verbose>0) {printf("%s()\n", __func__); fflush(stdout);}
 
  if(img==NULL || imgHasData(img)==0 || dsr==NULL) return(TPCERROR_FAIL);
 
  /* Check file format */
  switch(img->oformat) {
    case IMG_FORMAT_NIFTI_1D:
      dsr->n=1;
      break;
    case IMG_FORMAT_NIFTI_1S:
      dsr->n=1;
      break;
    case IMG_FORMAT_NIFTI_2D:
      dsr->n=2;
      break;
    case IMG_FORMAT_NIFTI_2S:
      dsr->n=2;
      break;
    default:
      return(TPCERROR_INVALID_FORMAT);
  }
 
  /* Find pixel value range */
  float pmin, pmax;
  int ret=imgMinMax(img, &pmin, &pmax);
  if(ret!=TPCERROR_OK) return(ret);
  if(verbose>1) {printf("  min=%g\n  max=%g\n", pmin, pmax); fflush(stdout);}
 
  /* Set NIfTI byte order to current machines byte order */
  dsr->byte_order=endianLittle();
 
  /* Initiate header structures with zeroes */
  memset(&dsr->h1, 0, sizeof(NIFTI_1_HEADER));
  memset(&dsr->h2, 0, sizeof(NIFTI_2_HEADER));
  memset(&dsr->e, 0, sizeof(NIFTI_EXTENDER));
 
 
  /* Set header */
  if(dsr->n==1) { // NIfTI-1
 
    dsr->h1.sizeof_hdr=NIFTI1_HEADER_SIZE; // 348
    strcpy(dsr->h1.data_type, ""); // not used in NIfTI
    strcpy(dsr->h1.db_name, ""); // not used in NIfTI
    dsr->h1.extents=16384; // not used in NIfTI, but required for Analyze compatibility
    dsr->h1.regular='r'; // not used in NIfTI, but required for Analyze compatibility
    dsr->h1.dim_info='\0'; // MRI slice ordering
    if(img->oformat==IMG_FORMAT_NIFTI_1S) {
      dsr->h1.vox_offset=NIFTI1_HEADER_SIZE+NIFTI1_HEADER_EXTENDER_SIZE; // 352
      strcpy(dsr->h1.magic, "n+1");
    } else {
      dsr->h1.vox_offset=0;
      strcpy(dsr->h1.magic, "ni1");
    }
    strcpy(dsr->h1.aux_file, "");
 
    /* Image dimension */
    for(int i=0; i<8; i++) dsr->h1.dim[i]=1;
    dsr->h1.dim[0]=4;
    dsr->h1.dim[1]=img->dimx;
    dsr->h1.dim[2]=img->dimy;
    dsr->h1.dim[3]=img->dimz;
    dsr->h1.dim[4]=img->dimt;
 
    /* Intent */
    dsr->h1.intent_code=NIFTI_INTENT_NONE;
    dsr->h1.intent_p1=0.0;
    dsr->h1.intent_p2=0.0;
    dsr->h1.intent_p3=0.0;
    strcpy(dsr->h1.intent_name, "");
 
    /* Save data as floats, so that there is no need to scale */
    dsr->h1.datatype=NIFTI_DT_FLOAT;
    dsr->h1.bitpix=32; // bits per pixel
    dsr->h1.scl_slope=1.0; // data as floats, so no need to scale
    dsr->h1.scl_inter=0.0; // data as floats, so no need to scale
    dsr->h1.glmax=pmax; // unused in NIfTI
    dsr->h1.glmin=pmin; // unused in NIfTI
    dsr->h1.cal_max=pmax;
    dsr->h1.cal_min=0.0; // scale display colours to have black at zero
 
    /* Pixel size */
    dsr->h1.slice_start=0;
    for(int i=0; i<8; i++) dsr->h1.pixdim[i]=0.0;
    // https://nifti.nimh.nih.gov/nifti-1/documentation/nifti1fields/nifti1fields_pages/qsform.html
    dsr->h1.pixdim[0]=1.0; // Set to either 1.0 or -1.0; default should be 1.0.
    dsr->h1.pixdim[1]=img->sizex;
    dsr->h1.pixdim[2]=img->sizey;
    dsr->h1.pixdim[3]=img->sizez;
    dsr->h1.slice_end=0;
    dsr->h1.slice_code=0;
    dsr->h1.xyzt_units=NIFTI_UNITS_MM; // Voxel size in mm in IMG structure
    dsr->h1.slice_duration=0.0;
    dsr->h1.toffset=0.0;
 
    /* Coordinate system */
    dsr->h1.qform_code=img->xform[0];
    dsr->h1.sform_code=img->xform[1];
    dsr->h1.quatern_b=img->quatern[0];
    dsr->h1.quatern_c=img->quatern[1];
    dsr->h1.quatern_d=img->quatern[2];
    dsr->h1.qoffset_x=img->quatern[3];
    dsr->h1.qoffset_y=img->quatern[4];
    dsr->h1.qoffset_z=img->quatern[5];
 
    /* Save pixel units into description (not standard NIfTI) */
    if(img->cunit==UNIT_UNKNOWN) strcpy(dsr->h1.descrip, "");
    else strlcpy(dsr->h1.descrip, unitName(img->cunit), 80);
 
    for(int i=0; i<4; i++) dsr->h1.srow_x[i]=img->srow[i];
    for(int i=0; i<4; i++) dsr->h1.srow_y[i]=img->srow[4+i];
    for(int i=0; i<4; i++) dsr->h1.srow_z[i]=img->srow[8+i];
 
    /* Extension is left as 0 0 0 0 */
 
  } else { // NIfTI-2
 
    /* Header size. Used for identifying the format and endianness of data */ 
    dsr->h2.sizeof_hdr=NIFTI2_HEADER_SIZE; // 540
 
    /* Magic string. Used for identifying the format and verifying validity */
    if(img->oformat==IMG_FORMAT_NIFTI_2S) strcpy(dsr->h2.magic, "n+2");
    else strcpy(dsr->h2.magic, "ni2");
    dsr->h2.magic[3]='\0';
    dsr->h2.magic[4]='\r';
    dsr->h2.magic[5]='\n';
    dsr->h2.magic[6]=(char)32;
    dsr->h2.magic[7]='\n';
 
    /* Save data as 32-bit floats, so that there is no need to scale */
    dsr->h2.datatype=NIFTI_DT_FLOAT;
    dsr->h2.bitpix=32; // bits per pixel
 
    /* Image dimensions */
    for(int i=0; i<8; i++) dsr->h2.dim[i]=1;
    dsr->h2.dim[0]=4;
    dsr->h2.dim[1]=img->dimx;
    dsr->h2.dim[2]=img->dimy;
    dsr->h2.dim[3]=img->dimz;
    dsr->h2.dim[4]=img->dimt;
 
    /* Intent */
    dsr->h2.intent_code=NIFTI_INTENT_NONE;
    dsr->h2.intent_p1=0.0;
    dsr->h2.intent_p2=0.0;
    dsr->h2.intent_p3=0.0;
    strcpy(dsr->h2.intent_name, "");
 
    /* Grid spacings */
    dsr->h2.pixdim[0]=1.0; // Set to either 1.0 or -1.0; default should be 1.0.
    dsr->h2.pixdim[1]=img->sizex;
    dsr->h2.pixdim[2]=img->sizey;
    dsr->h2.pixdim[3]=img->sizez;
    for(int i=4; i<8; i++) dsr->h2.pixdim[i]=0.0;
 
    /* Offset into pixel data; zero for dual file */
    if(img->oformat==IMG_FORMAT_NIFTI_2S) {
      dsr->h2.vox_offset=NIFTI2_HEADER_SIZE+NIFTI2_HEADER_EXTENDER_SIZE; // 544
    } else {
      dsr->h2.vox_offset=0; 
    }
 
    /* Scaling */
    dsr->h2.scl_slope=1.0; // data as floats, so no need to scale
    dsr->h2.scl_inter=0.0; // data as floats, so no need to scale
    dsr->h2.cal_max=pmax;
    dsr->h2.cal_min=0.0; // scale display colours to have black at zero
 
    dsr->h2.slice_duration=0.0;
    dsr->h2.toffset=0.0;
    dsr->h2.slice_start=0;
    dsr->h2.slice_end=img->dimz-1;
 
    /* Save pixel units into description (not standard NIfTI) */
    if(img->cunit==UNIT_UNKNOWN) strcpy(dsr->h2.descrip, "");
    else strlcpy(dsr->h2.descrip, unitName(img->cunit), 80);
 
    strcpy(dsr->h2.aux_file, "");
 
    /* Coordinate system */
    dsr->h2.qform_code=img->xform[0];
    dsr->h2.sform_code=img->xform[1];
    dsr->h2.quatern_b=img->quatern[0];
    dsr->h2.quatern_c=img->quatern[1];
    dsr->h2.quatern_d=img->quatern[2];
    dsr->h2.qoffset_x=img->quatern[3];
    dsr->h2.qoffset_y=img->quatern[4];
    dsr->h2.qoffset_z=img->quatern[5];
    for(int i=0; i<4; i++) dsr->h2.srow_x[i]=img->srow[i];
    for(int i=0; i<4; i++) dsr->h2.srow_y[i]=img->srow[4+i];
    for(int i=0; i<4; i++) dsr->h2.srow_z[i]=img->srow[8+i];
 
    dsr->h2.slice_code=0;
    dsr->h2.xyzt_units=NIFTI_UNITS_MM+NIFTI_UNITS_SEC;
    dsr->h2.dim_info='\0'; // MRI slice ordering
 
    /* Extension is left as 0 0 0 0 */
 
  }
 
  return(TPCERROR_OK);
}

Referenced by imgWriteNifti().

imgToSIF()

int imgToSIF ( IMG * img, TAC * sif, int copy_header, int copy_frames, int copy_counts, int verbose )

extern

Set SIF contents based on data in IMG.

Todo

Add tests.

See also

imgFromSIF, imgInit, tacInit, imgWrite, tacWrite, tacWriteSIF

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to IMG struct from which content is copied to SIF.
sif	Pointer to SIF, which actually is stored as a TAC format. Previous contents are optionally kept.

Precondition

Must be initialized with tacInit().

Parameters

copy_header	Select whether header contents are copied (1) or not copied (0) to SIF.
copy_frames	Select whether frame times are copied (1) or not copied (0) to SIF.
copy_counts	Select whether counts are copied (1) or not copied (0) to SIF, or created if they do not exist (2).
verbose	Verbose level; if zero, then nothing is printed to stderr or stdout.

Definition at line 18 of file imagesif.c.

  {
  if(verbose>0) {
    printf("%s(img, sif, %d, %d, %d, ...)\n", __func__, copy_header, copy_frames, copy_counts);
    fflush(stdout);
  }
 
  if(img==NULL || sif==NULL) return(TPCERROR_FAIL);
 
  /* Verify that IMG contains frame times */
  if(!imgHasTimes(img)) {
    if(verbose>0) printf("  image does not contain frame times.\n");
    /* If not, then frame times cannot be copied */
    copy_frames=0;
    /* and counts can not be created */
    if(copy_counts==2) copy_counts=1;
  }
 
  /* Check if count data needs to be created */
  if(copy_counts==2 && imgHasCounts(img)!=0) copy_counts=1;
 
  /* Verify that IMG contains isotope information */
  if(img->isot==ISOTOPE_UNKNOWN) {
    if(verbose>0) printf("  image does not contain isotope halflife.\n");
    /* not, then count data can not be created */
    if(copy_counts==2) copy_counts=1;
  }
 
  /* Allocate memory for SIF if necessary */
  int ret=0;
  if(sif->sampleNr<1 || sif->tacNr<1) {
    ret=tacAllocate(sif, img->dimt, 2);
    sif->sampleNr=img->dimt; sif->tacNr=2;
  } else if((copy_frames || copy_counts) && sif->sampleNr!=img->dimt) {
    if(sif->sampleNr<img->dimt)
      ret=tacAllocateMoreSamples(sif, img->dimt-sif->sampleNr);
    sif->sampleNr=img->dimt;
  }
  if((copy_frames || copy_counts) && sif->tacNr<2) {
    ret=tacAllocateMore(sif, 2-sif->tacNr);
    sif->tacNr=2;
  }
  if(ret!=0) return(TPCERROR_OUT_OF_MEMORY);
 
  /* Copy SIF header */
  if(copy_header) {
    if(verbose>1) printf("  copying header fields.\n");
    tacSetHeaderScanstarttime(&sif->h, img->scanStart);
    tacSetHeaderStudynr(&sif->h, img->studyNr);
    tacSetIsotope(sif, img->isot);
  }
 
  /* Copy frame times */
  if(copy_frames) {
    if(verbose>1) printf("  copying frame times.\n");
    for(int i=0; i<img->dimt; i++) {
      sif->x1[i]=img->x1[i]; sif->x2[i]=img->x2[i]; sif->x[i]=img->x[i];
    }
    sif->tunit=UNIT_SEC;
  }
 
  /* Copy counts */
  if(copy_counts==1) {
    if(verbose>1) printf("  copying count data.\n");
    for(int i=0; i<img->dimt; i++) {
      sif->c[0].y[i]=img->prompts[i];
      sif->c[1].y[i]=img->randoms[i];
    }
  }
 
  /* Create counts, if required and possible */
  if(copy_counts==2) {
    if(verbose>1) printf("  creating count data.\n");
 
    fprintf(stderr, "Warning: imgToSIF() does not yet guess counts.\n"); fflush(stderr);
 
    for(int i=0; i<img->dimt; i++) {
      sif->c[0].y[i]=0.0;
      sif->c[1].y[i]=0.0;
    }
  }
 
  return(TPCERROR_OK);
}

Referenced by imgWriteNifti().

imgWrite()

int imgWrite ( IMG * img, const char * fname, TPCSTATUS * status )

extern

Write image file from IMG structure.

Remarks

Stub function.

See also

imgInit, imgRead, imgFree

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to image structure.
fname	Pointer to the file name, possibly path in case of DICOM.
status	Pointer to status data; enter NULL if not needed.

Definition at line 214 of file imageio.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(img, %s)\n", __func__, fname); fflush(stdout);}
 
  if(img==NULL || strnlen(fname, 2)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
  if(!imgHasData(img)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(TPCERROR_NO_DATA);
  }
 
  /* Based on file name, guess file format, and use it, or check that specified format is ok */
  imgformat format=imgFormatFromFName(fname, NULL);
  if(format!=IMG_FORMAT_UNKNOWN) {
    if(img->oformat==IMG_FORMAT_UNKNOWN) {
      img->oformat=format;
      if(verbose>1)
        printf("based on file name, output format set to %s\n", imgFormatDescr(img->oformat));
    } else if(format!=img->oformat) {
      if(verbose>1)
        printf("file name is not suitable for format %s\n", imgFormatDescr(img->oformat));
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
      return(TPCERROR_INVALID_FILENAME);
    }
  }
  /* If output format is not yet set, then use input format */
  if(img->oformat==IMG_FORMAT_UNKNOWN && img->format!=IMG_FORMAT_UNKNOWN) {
    img->oformat=img->format;
    if(verbose>1) printf("output set to same as input format (%s)\n", imgFormatDescr(img->oformat));
  }
  /* If output format still unknown, set it to singe-file NIfTI */
  if(img->oformat==IMG_FORMAT_UNKNOWN) {
    img->oformat=IMG_FORMAT_NIFTI_1S;
    if(verbose>1) printf("output format set to default (%s)\n", imgFormatDescr(img->oformat));
  }
 
  /* Fill output header */
  if(imgFillOHeader(img, status)!=TPCERROR_OK) return(TPCERROR_INVALID_HEADER);
 
  /* Write image */
  if(img->oformat==IMG_FORMAT_DICOM) {
 
    if(verbose>1) printf("writing %s\n", imgFormatDescr(img->oformat));
 
  } else if(img->oformat==IMG_FORMAT_NIFTI_1D || img->oformat==IMG_FORMAT_NIFTI_1S) {
 
    if(verbose>1) printf("writing %s\n", imgFormatDescr(img->oformat));
    return(imgWriteNifti(img, fname, status));
 
  } else if(img->oformat==IMG_FORMAT_E7) {
 
    if(verbose>1) printf("writing %s\n", imgFormatDescr(img->oformat));
 
    /* Create main header data block */
    unsigned char mh[ECATBLKSIZE];
    int ret=ecatWriteMainheader(&img->oh, mh, status);
    if(ret!=TPCERROR_OK) return(ret);
 
    /* Open file, write main header block, and initiate matrix list */
 
    /* Process each frame */
    for(unsigned int ti=0; ti<img->dimt; ti++) {
      /* Create subheader data block */
 
      /* Create blocks for pixel data */
 
      /* Write header and data blocks */
 
    }
 
    /* Close file */
 
  } else if(img->oformat==IMG_FORMAT_E63) {
 
    if(verbose>1) printf("writing %s\n", imgFormatDescr(img->oformat));
 
    /* Create main header data block */
 
    /* Open file, write main header block, and initiate matrix list */
 
    /* Process each plane and frame */
    for(unsigned int ti=0; ti<img->dimt; ti++) {
      for(unsigned int zi=0; zi<img->dimz; zi++) {
        /* Create subheader data block */
 
        /* Create blocks for pixel data */
 
        /* Write header and data blocks */
 
      }
    }
 
    /* Close file */
 
  } else if(img->oformat==IMG_FORMAT_FLAT) {
 
    if(verbose>1) printf("writing %s\n", imgFormatDescr(img->oformat));
    FILE *fp=NULL;
    float *fdata=NULL, *fptr;
    if(verbose>2) printf("opening %s\n", fname);
    if((fp=fopen(fname, "wb")) == NULL) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
      return(TPCERROR_CANNOT_WRITE);
    }
    /* Allocate memory for float data for one frame */
    unsigned long long pxlNr=img->dimz*img->dimy*img->dimx;
    if(verbose>2) printf("allocating memory for %llu pixels\n", pxlNr);
    fdata=(float*)malloc(pxlNr*sizeof(float));
    if(fdata==NULL) {
      fclose(fp);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
      return(TPCERROR_OUT_OF_MEMORY);
    }
    /* Write one frame at a time */
    if(verbose>2) printf("writing frame at a time\n");
    for(int ti=0; ti<img->dimt; ti++) {
      fptr=fdata;
      for(int zi=0; zi<img->dimz; zi++)
        for(int yi=0; yi<img->dimy; yi++)
          for(int xi=0; xi<img->dimx; xi++)
            *fptr++=img->m[zi][yi][xi][ti];
      if(fwrite((float*)fdata, 4, pxlNr, fp) != pxlNr) {
        fclose(fp); remove(fname); free(fdata);
        statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
        return(TPCERROR_CANNOT_WRITE);
      }
    } /* next frame */
    /* Close file, free memory */
    fclose(fp); free(fdata);
    if(verbose>2) printf("binary file written\n");
    /* Write matrix information file; no error if not successful. */
    char miffile[strlen(fname)+5]; sprintf(miffile, "%s.mif", fname);
    if(verbose>2) printf("opening matrix information file %s\n", miffile);
    if((fp=fopen(miffile, "w")) == NULL) {
      if(verbose>0) fprintf(stderr, "Error: cannot open MIF for write.\n");
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
      return(TPCERROR_OK);
    }
    if(verbose>2) printf("writing matrix information file %s\n", miffile);
    int n=fprintf(fp, "%d %d %d %d\n", img->dimz, img->dimt, img->dimx, img->dimy);
    fclose(fp);
    if(n<7) {
      if(verbose>0) fprintf(stderr, "Error: cannot write in file %s\n", miffile); 
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
      return(TPCERROR_OK);
    }
    if(verbose>2) printf("matrix information saved in %s\n", miffile);
 
  } else {
    if(verbose>0) printf("output format (%s) not supported\n", imgFormatDescr(img->oformat));
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_UNSUPPORTED);
    return(TPCERROR_UNSUPPORTED);
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

imgWriteDICOM()

int imgWriteDICOM ( IMG * img, const char * fname, TPCSTATUS * status )

extern

Write DICOM image into IMG data structure.

Remarks

Stub function.

Todo

Add tests.

See also

imgInit, imgRead, imgWrite, imgFree, imgReadDICOM

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to image structure.
fname	Pointer to the name of DICOM folder. Must be empty or not exist.
status	Pointer to status data; obligatory.

Definition at line 335 of file imagedcm.c.

  {
  if(status==NULL) return(TPCERROR_FAIL);
  if(status->verbose>1) {printf("%s(..., %s, ...)\n", __func__, fname); fflush(stdout);}
  if(img==NULL || img->dimt<1 || img->dimz<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA); return(status->error);}
  if(fname==NULL || strnlen(fname, 1)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME); return(status->error);}
 
  /*
   *  Make DCMFILE data structure
   */
  DCMFILE dcm; dcmfileInit(&dcm);
  /* Set Transfer UID */
  dcm.truid=DCM_TRUID_LEE;
  {
    DCMTAG tag;
    unsigned int ul;
    char buf[256];
 
    /* Meta Element Group */
    dcmTagSet(&tag, 0x0002, 0x0000); ul=0;
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UL, 0xFFFFFFFF, (char*)&ul, 0);
    dcmTagSet(&tag, 0x0002, 0x0001); // File Meta Information Version
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_OB, 2, "\0\1", 0);
    dcmTagSet(&tag, 0x0002, 0x0002); // Media Stored SOP Class UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.2.840.10008.5.1.4.1.1.128", 0);
    dcmTagSet(&tag, 0x0002, 0x0003); // Media Stored SOP Instance UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.2.826.0.1.3417726.3.0.20180616", 0);
    dcmTagSet(&tag, 0x0002, 0x0010); // Transfer Syntax UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, dcmTrUIDString(dcm.truid), 0);
    dcmTagSet(&tag, 0x0002, 0x0012); // Implementation Class UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.3.6.1.4.1.25403.1.1.1", 0);
    dcmTagSet(&tag, 0x0002, 0x0013); // Implementation Version Name
    snprintf(buf, 256, "TPCCLIB_%d.%d.%d", tpcclib_VERSION_MAJOR, tpcclib_VERSION_MINOR, tpcclib_VERSION_PATCH);
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_SH, 0xFFFFFFFF, buf, 0);
 
    /* Identifying Group */
    dcmTagSet(&tag, 0x0008, 0x0000); ul=0;
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UL, 0xFFFFFFFF, (char*)&ul, 0);
    dcmTagSet(&tag, 0x0008, 0x0008); // Image Type
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_CS, 17, "ORIGINAL/PRIMARY", 0);
    {
      time_t t=time(NULL);
      struct tm *local_tm; local_tm=localtime(&t);
      dcmTagSet(&tag, 0x0008, 0x0012); // Instance Creation Date
      strftime(buf, 256, "%Y%m%d", local_tm);
      dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_DA, 0xFFFFFFFF, buf, 0);
      dcmTagSet(&tag, 0x0008, 0x0013); // Instance Creation Time
      strftime(buf, 256, "%H%M%S", local_tm);
      dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_DA, 0xFFFFFFFF, buf, 0);
    }
    dcmTagSet(&tag, 0x0008, 0x0016); // SOP Class UID; see dcmsop.c
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.2.840.10008.5.1.4.1.1.128", 0);
    dcmTagSet(&tag, 0x0008, 0x0018); // SOP Instance UID
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UI, 0xFFFFFFFF, "1.2.826.0.1.3417726.3.0.20180616", 0);
    dcmTagSet(&tag, 0x0008, 0x0060); // Modality
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_CS, 3, "PT", 0);
 
    /* Acquisition Group */
    dcmTagSet(&tag, 0x0018, 0x0000); ul=0;
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UL, 0xFFFFFFFF, (char*)&ul, 0);
    dcmTagSet(&tag, 0x0018, 0x0050);  // Slice thickness (mm)
    snprintf(buf, 256, "%g", img->sizez);
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_DS, 0xFFFFFFFF, buf, 0);
    dcmTagSet(&tag, 0x0018, 0x1075); // half-life
    if(img->isot!=ISOTOPE_UNKNOWN) {
      snprintf(buf, 256, "%g", 60.0*isotopeHalflife(img->isot));
      dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_DS, 0xFFFFFFFF, buf, 0);
    }
 
    /* Relationship Group */
    dcmTagSet(&tag, 0x0020, 0x0000); ul=0;
    dcmAddItem(&dcm, NULL, 0, tag, DCM_VR_UL, 0xFFFFFFFF, (char*)&ul, 0);
 
 
 
  }
 
 
  /* Create folder, or if it already exists check that it is empty */
  if(pathExist(fname)) {
    if(status->verbose>2) {
      fprintf(stdout, "  checking that folder '%s' is empty\n", fname); fflush(stdout);}
    if(pathFileNr(fname)>0) {
      dcmfileFree(&dcm);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FILE_EXISTS);
      return(status->error);
    }
  } else {
    if(status->verbose>2) {
      fprintf(stdout, "  creating subdirectory '%s'\n", fname); fflush(stdout);}
    if(pathCreate(fname)) {
      dcmfileFree(&dcm);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
      return(status->error);
    }
  }
 
  /*
   *  Write each image plane in its own DICOM file
   */
  for(int zi=0; zi<img->dimz; zi++) {
    char dname[strlen(fname)+25];
    snprintf(dname, strlen(fname)+25, "%s/%0*d.dcm", fname, 1+(int)floor(log10(img->dimz)), 1+zi);
    if(status->verbose>2) {printf("  writing '%s'\n", dname); fflush(stdout);}
 
    /* Write */
    dcmFileWrite(dname, &dcm, status);
    if(status->error!=TPCERROR_OK) {
      dcmfileFree(&dcm);
      dcmFileRemove(dname, NULL);
      pathRemove(fname);
      return(status->error);
    }
  } // next plane
  dcmfileFree(&dcm);
 
  if(0) { // Remove files and empty folder in case of error
    dcmFileRemove(fname, NULL);
    pathRemove(fname);
  }
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

imgWriteNifti()

int imgWriteNifti ( IMG * img, const char * fname, TPCSTATUS * status )

extern

Write NIfTI from IMG data structure.

IMG field oformat (or secondarily format) determines whether NIfTI is written in single file format (*.nii) or dual file format (*.hdr and *.img).

SIF file is saved to store frame times, if frame times are available.

Todo

Add tests.

See also

imgInit, imgRead, imgWrite, imgFree, niftiCreateFNames, niftiWriteHeader

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to image structure. Any previous contents removed.
fname	Pointer to the database name. Path not accepted.
status	Pointer to status data; enter NULL if not needed.

Definition at line 390 of file imagenii.c.

  {
  int verbose=0; if(status!=NULL) verbose=status->verbose;
  if(verbose>0) {printf("%s(..., %s, ...)\n", __func__, fname); fflush(stdout);}
 
  if(img==NULL || strnlen(fname, 2)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_FAIL);
    return(TPCERROR_FAIL);
  }
  if(!imgHasData(img)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_NO_DATA);
    return(TPCERROR_NO_DATA);
  }
 
  /* Determine the output format */
  if(verbose>1) {printf("determine output format\n"); fflush(stdout);}
  if(img->oformat==IMG_FORMAT_UNKNOWN) {
    char *cptr=filenameGetExtension(fname);
    if(cptr!=NULL && strcasecmp(cptr, ".nii")==0) img->oformat=IMG_FORMAT_NIFTI_1S;
    else img->oformat=img->format;
  }
  if(img->oformat!=IMG_FORMAT_NIFTI_1D && img->oformat!=IMG_FORMAT_NIFTI_1S &&
     img->oformat!=IMG_FORMAT_NIFTI_2D && img->oformat!=IMG_FORMAT_NIFTI_2S) 
    img->oformat=IMG_FORMAT_NIFTI_1S;
 
  /* Get database name, in case file name was given with extension */
  if(verbose>1) {printf("get dbname\n"); fflush(stdout);}
  char basename[FILENAME_MAX];
  strlcpy(basename, fname, FILENAME_MAX); niftiBasename(basename);
  /* Check that there is something left after path */
  char basenamewopath[FILENAME_MAX];
  strcpy(basenamewopath, basename); filenameRmPath(basenamewopath);
  if(strlen(basenamewopath)<1) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return(TPCERROR_INVALID_FILENAME);
  }
  /* Create output file names (now to catch some errors) */
  if(verbose>1) {printf("determine output file names\n"); fflush(stdout);}
  char hdrfile[FILENAME_MAX], imgfile[FILENAME_MAX], siffile[FILENAME_MAX];
  hdrfile[0]=imgfile[0]=siffile[0]=(char)0;
  if(niftiCreateFNames(basename, hdrfile, imgfile, siffile, img->oformat)!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return(TPCERROR_INVALID_FILENAME);
  }
  /* If we have path, then create it */
  char filepath[FILENAME_MAX];
  strcpy(filepath, basename); filenameRmFile(filepath);
  if(pathCreate(filepath)) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_FILENAME);
    return(TPCERROR_INVALID_FILENAME);
  }
 
 
  /* Fill the header structure */
  NIFTI_DSR dsr;
  int ret=imgSetNiftiHeader(img, &dsr, verbose-1);
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
 
  /* Delete previous NIfTI */
  if(verbose>1) {printf("removing any previous files\n"); fflush(stdout);}
  if(fileExist(hdrfile)) remove(hdrfile);
  if(fileExist(imgfile)) remove(imgfile);
  //if(fileExist(siffile)) remove(siffile); // NO, not this!
 
  /* Write the header */
  if(verbose>1) {printf("writing the header\n"); fflush(stdout);}
  if(img->oformat==IMG_FORMAT_NIFTI_1D || img->oformat==IMG_FORMAT_NIFTI_2D) {
    ret=niftiWriteHeader(hdrfile, &dsr, verbose-1);
  } else if(img->oformat==IMG_FORMAT_NIFTI_1S || img->oformat==IMG_FORMAT_NIFTI_2S) {
    ret=niftiWriteHeader(imgfile, &dsr, verbose-1);
  }
  if(ret!=TPCERROR_OK) {
    statusSet(status, __func__, __FILE__, __LINE__, ret);
    return(ret);
  }
  /* Write the pixel data */
  if(verbose>1) {printf("opening file for pixel data\n"); fflush(stdout);}
  FILE *fp;
  if(img->oformat==IMG_FORMAT_NIFTI_1D || img->oformat==IMG_FORMAT_NIFTI_2D) {
    fp=fopen(imgfile, "wb");
  } else { // if(img->oformat==IMG_FORMAT_NIFTI_1S || if(img->oformat==IMG_FORMAT_NIFTI_2S
    fp=fopen(imgfile, "r+b");
  }
  if(fp==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_OPEN);
    return(TPCERROR_CANNOT_OPEN);
  }
  /* Move pointer to the start of pixel data */
  long int pos;
  if(dsr.n==1) pos=(int)dsr.h1.vox_offset; else pos=(int)dsr.h2.vox_offset;
  if(fseek(fp, pos, SEEK_SET)!=0) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    fclose(fp); return(TPCERROR_CANNOT_WRITE);
  }
  /* Allocate memory for pixel matrix */
  unsigned long long voxNr=img->dimt*img->dimz*img->dimy*img->dimx;
  float *fdata;
  fdata=(float*)calloc(voxNr, sizeof(float));
  if(fdata==NULL) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OUT_OF_MEMORY);
    fclose(fp); return(TPCERROR_OUT_OF_MEMORY);
  }
  /* Write voxel values as floats */
  float *fptr=fdata;
  for(unsigned int fi=0; fi<img->dimt; fi++)
    for(unsigned int zi=0; zi<img->dimz; zi++)
      for(unsigned int yi=0; yi<img->dimy; yi++)
        for(unsigned int xi=0; xi<img->dimx; xi++, fptr++)
          *fptr=img->m[zi][yi][xi][fi];
  if(verbose>1) {printf("writing pixel data\n"); fflush(stdout);}
  fptr=fdata;
  if(fwrite(fptr, sizeof(float), voxNr, fp) != voxNr) {
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
    free(fdata); fclose(fp); return(TPCERROR_CANNOT_WRITE);
  }
 
  free(fdata);
  fclose(fp);
 
  /* Frame times into SIF */
  if(!imgHasTimes(img)) { // no frame times to save
    if(verbose>1) printf("  no frame times to save into SIF\n");
    statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
    return(TPCERROR_OK);
  }
  {
    if(verbose>1) {printf("making SIF\n"); fflush(stdout);}
    TAC sif; tacInit(&sif);
    /* Try to read existing SIF */
    if(fileExist(siffile) && tacRead(&sif, siffile, status)==TPCERROR_OK && sif.sampleNr==img->dimt)
    {
      /* If SIF was found and frameNr matches with image, then update the SIF contents, 
         but keep counts, in case previous SIF comes with actual count info from scanner */
        ret=imgToSIF(img, &sif, 1, 1, 0, verbose-3);
    } else {
      /* otherwise create SIF contents */
      ret=imgToSIF(img, &sif, 1, 1, 2, verbose-3);
    }
    if(ret!=TPCERROR_OK) {
      if(verbose>0) fprintf(stderr, "  Error: cannot create SIF contents.\n");
      tacFree(&sif);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_INVALID_X);
      return(TPCERROR_INVALID_X);
    }
    if(verbose>1) {printf("writing SIF\n"); fflush(stdout);}
    fp=fopen(siffile, "w");
    if(fp==NULL) {
      tacFree(&sif);
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_OPEN);
      return(TPCERROR_CANNOT_OPEN);
    }
    ret=tacWrite(&sif, fp, TAC_FORMAT_SIF, 0, status);
    fclose(fp); tacFree(&sif);
    if(ret!=TPCERROR_OK) {
      statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_CANNOT_WRITE);
      return(TPCERROR_CANNOT_WRITE);
    }
  }
 
 
  statusSet(status, __func__, __FILE__, __LINE__, TPCERROR_OK);
  return(TPCERROR_OK);
}

Referenced by imgWrite().

imgXRange()

int imgXRange ( IMG * img, double * xmin, double * xmax )

extern

Get the range of x values (times) in the image data.

Data is not modified. Data does not need to be sorted.

See also

imgMinMax, imgHasData

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to IMG structure.
xmin	Pointer to variable for min x value (NULL if not needed).
xmax	Pointer to variable for max x value (NULL if not needed).

Definition at line 430 of file image.c.

  {
  if(xmin!=NULL) *xmin=nan("");
  if(xmax!=NULL) *xmax=nan("");
  if(img==NULL || !imgHasData(img)) return(TPCERROR_FAIL);
  float f1=nanf(""), f2=nanf("");
  for(int i=0; i<img->dimt; i++) {
    if(isfinite(img->x1[i])) {
      if(isnan(f1) || img->x1[i]<f1) f1=img->x1[i];
    }
    if(isfinite(img->x2[i])) {
      if(isnan(f2) || img->x2[i]>f2) f2=img->x2[i];
    }
  }
  if(xmin!=NULL) *xmin=f1;
  if(xmax!=NULL) *xmax=f2;
  if(isnan(f1) || isnan(f2)) return(TPCERROR_FAIL);
  return(TPCERROR_OK);
}

imgXUnitConvert()

int imgXUnitConvert	(	IMG *	img,
		const int	u )

Convert X values (frame times) in IMG structure to the specified unit.

See also

tacXUnitConvert, unitConversionFactor

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

img	Pointer to IMG structure.
u	Enum unit

Definition at line 462 of file image.c.

  {
  if(img==NULL || img->dimt<1) return(TPCERROR_NO_DATA);
 
  /* If both units are unknown, then return suitable error code */
  if(u==UNIT_UNKNOWN && img->tunit==UNIT_UNKNOWN) return(TPCERROR_UNKNOWN_UNIT);
 
  /* Determine the correction factor */
  double cf;
  cf=unitConversionFactor(img->tunit, u);
  if(isnan(cf)) return(TPCERROR_INCOMPATIBLE_UNIT);
  
  /* Multiply x values */
  for(int i=0; i<img->dimt; i++) {
    img->x[i]*=cf; img->x1[i]*=cf; img->x2[i]*=cf;
  }
  /* Set new unit in the header */
  img->tunit=u;
 
  return(TPCERROR_OK);
}

niftiCreateFNames()

int niftiCreateFNames ( const char * filename, char * hdrfile, char * imgfile, char * siffile, int fileformat )

extern

Construct the file names for NIfTI image.

See also

imgReadNifti, imgWriteNifti, niftiRemoveFNameExtension

Returns

enum tpcerror (TPCERROR_OK when successful).

Parameters

filename	Filename, either header file, image file, or base name without extensions, but possibly with path name. This string is never modified.
hdrfile	Header filename will be written in this char pointer (space needs to allocated by caller); in single file format this will be set to the name of the single file; enter NULL if not needed.
imgfile	Image filename will be written in this char pointer (space needs to allocated by caller); in single file format this will be set to the name of the single file; enter NULL if not needed.
siffile	SIF filename will be written in this char pointer (space needs to allocated by caller); enter NULL if not needed.
fileformat	NIfTI file format, either IMG_FORMAT_NIFTI_1D, IMG_FORMAT_NIFTI_1S, IMG_NIFTI_FORMAT_2D, or IMG_FORMAT_NIFTI_2S.

Definition at line 17 of file imagenii.c.

  {
  if(hdrfile!=NULL) strcpy(hdrfile, "");
  if(imgfile!=NULL) strcpy(imgfile, "");
  if(siffile!=NULL) strcpy(siffile, "");
  if(filename==NULL) return(TPCERROR_INVALID_FILENAME);
 
  char basename[FILENAME_MAX];
  strlcpy(basename, filename, FILENAME_MAX); niftiBasename(basename);
  /* Check that there is something left after path */
  char basenamewopath[FILENAME_MAX];
  strcpy(basenamewopath, basename); filenameRmPath(basenamewopath);
  if(strlen(basenamewopath)<1) return(TPCERROR_INVALID_FILENAME);
 
  /* Create database filenames */
  if(fileformat==IMG_FORMAT_NIFTI_1D || fileformat==IMG_FORMAT_NIFTI_2D) {
    if(hdrfile!=NULL) snprintf(hdrfile, FILENAME_MAX, "%s.hdr", basename);
    if(imgfile!=NULL) snprintf(imgfile, FILENAME_MAX, "%s.img", basename);
  } else if(fileformat==IMG_FORMAT_NIFTI_1S || fileformat==IMG_FORMAT_NIFTI_2S) {
    if(hdrfile!=NULL) snprintf(hdrfile, FILENAME_MAX, "%s.nii", basename);
    if(imgfile!=NULL) snprintf(imgfile, FILENAME_MAX, "%s.nii", basename);
  } else {
    return(TPCERROR_INVALID_FORMAT);
  }
  if(siffile!=NULL) snprintf(siffile, FILENAME_MAX, "%s.sif", basename);
 
  return(TPCERROR_OK);
}

Referenced by imgWriteNifti().