nifti1_io.c

来自「DTMK软件开发包,此为开源软件,是一款很好的医学图像开发资源.」· C语言 代码 · 共 1,582 行 · 第 1/5 页

      nim->nw = nim->dim[7];
   nim->dw = nim->pixdim[7];

   for( c = 1, nim->nvox = 1; c <= nim->dim[0]; c++ )
      nim->nvox *= nim->dim[c];

   /* compute ndim, assuming it can be no larger than the old one */
   for( ndim = nim->dim[0]; (ndim > 1) && (nim->dim[ndim] <= 1); ndim-- )
       ;

   if( g_opts.debug > 2 ){
      fprintf(stderr,"+d ndim = %d -> %d\n",nim->ndim, ndim);
      fprintf(stderr," --> (%d,%d,%d,%d,%d,%d,%d)\n",
              nim->nx, nim->ny, nim->nz, nim->nt, nim->nu, nim->nv, nim->nw);
   }

   nim->dim[0] = nim->ndim = ndim;

   return 0;
}


/*----------------------------------------------------------------------*/
/*! Load the image data from disk into an already-prepared image struct.
 *
 * \param    nim      - initialized nifti_image, without data
 * \param    nbricks  - the length of blist (must be 0 if blist is NULL)
 * \param    blist    - an array of xyz volume indices to read (can be NULL)
 * \param    NBL      - pointer to struct where resulting data will be stored
 *
 * If blist is NULL, read all sub-bricks.
 * 
 * \return the number of loaded bricks (NBL->nbricks),
 *    0 on failure, < 0 on error
 *
 * NOTE: it is likely that another function will copy the data pointers
 *       out of NBL, in which case the only pointer the calling function
 *       will want to free is NBL->bricks (not each NBL->bricks[i]).
*//*--------------------------------------------------------------------*/
int nifti_image_load_bricks( nifti_image * nim , int nbricks,
                             const int * blist, nifti_brick_list * NBL )
{
   int     * slist = NULL, * sindex = NULL, rv;
   znzFile   fp;

   /* we can have blist == NULL */
   if( !nim || !NBL ){
      fprintf(stderr,"** nifti_image_load_bricks, bad params (%p,%p)\n",
              (void *)nim, (void *)NBL);
      return -1;
   }

   if( blist && nbricks <= 0 ){
      if( g_opts.debug > 1 )
         fprintf(stderr,"-d load_bricks: received blist with nbricks = %d,"
                        "ignoring blist\n", nbricks);
      blist = NULL; /* pretend nothing was passed */
   }

   if( blist && ! valid_nifti_brick_list(nim, nbricks, blist, g_opts.debug>0) )
      return -1;

   /* for efficiency, let's read the file in order */
   if( blist && nifti_copynsort( nbricks, blist, &slist, &sindex ) != 0 )
      return -1;

   /* open the file and position the FILE pointer */
   fp = nifti_image_load_prep( nim );
   if( !fp ){
      if( g_opts.debug > 0 )
         fprintf(stderr,"** nifti_image_load_bricks, failed load_prep\n");
      if( blist ){ free(slist); free(sindex); }
      return -1;
   }

   /* this will flag to allocate defaults */
   if( !blist ) nbricks = 0;
   if( nifti_alloc_NBL_mem( nim, nbricks, NBL ) != 0 ){
      if( blist ){ free(slist); free(sindex); }
      znzclose(fp);
      return -1;
   }

   rv = nifti_load_NBL_bricks(nim, slist, sindex, NBL, fp);

   if( rv != 0 ){
      nifti_free_NBL( NBL );  /* failure! */
      NBL->nbricks = 0; /* repetative, but clear */
   }

   if( slist ){ free(slist); free(sindex); }

   znzclose(fp);

   return NBL->nbricks;
}


/*----------------------------------------------------------------------*/
/*! nifti_free_NBL      - free all pointers and clear structure
 *
 * note: this does not presume to free the structure pointer
*//*--------------------------------------------------------------------*/
void nifti_free_NBL( nifti_brick_list * NBL )
{
   int c;

   if( NBL->bricks ){
      for( c = 0; c < NBL->nbricks; c++ )
         if( NBL->bricks[c] ) free(NBL->bricks[c]);
      free(NBL->bricks);
      NBL->bricks = NULL;
   }

   NBL->bsize = NBL->nbricks = 0;
}


/*----------------------------------------------------------------------
 * nifti_load_NBL_bricks      - read the file data into the NBL struct
 *
 * return 0 on success, -1 on failure
 *----------------------------------------------------------------------*/
static int nifti_load_NBL_bricks( nifti_image * nim , int * slist, int * sindex,
                                  nifti_brick_list * NBL, znzFile fp )
{
   size_t oposn, fposn;      /* orig and current file positions */
   size_t rv;
   long   test;
   int    c;
   int    prev, isrc, idest; /* previous and current sub-brick, and new index */

   test = znztell(fp);  /* store current file position */
   if( test < 0 ){
      fprintf(stderr,"** load bricks: ztell failed??\n");
      return -1;
   }
   fposn = oposn = test;

   /* first, handle the default case, no passed blist */
   if( !slist ){
      for( c = 0; c < NBL->nbricks; c++ ) {
         rv = nifti_read_buffer(fp, NBL->bricks[c], NBL->bsize, nim);
         if( rv != NBL->bsize ){
            fprintf(stderr,"** load bricks: cannot read brick %d from '%s'\n",
                    c, nim->iname ? nim->iname : nim->fname);
            return -1;
         }
      }
      if( g_opts.debug > 1 )
         fprintf(stderr,"+d read %d default %u-byte bricks from file %s\n",
                 NBL->nbricks, (unsigned int)NBL->bsize,
                 nim->iname ? nim->iname:nim->fname );
      return 0;
   }

   if( !sindex ){
      fprintf(stderr,"** load_NBL_bricks: missing index list\n");
      return -1;
   }

   prev = -1;   /* use prev for previous sub-brick */
   for( c = 0; c < NBL->nbricks; c++ ){
       isrc = slist[c];   /* this is original brick index (c is new one) */
       idest = sindex[c]; /* this is the destination index for this data */

       /* if this sub-brick is not the previous, we must read from disk */
       if( isrc != prev ){

          /* if we are not looking at the correct sub-brick, scan forward */
          if( fposn != (oposn + isrc*NBL->bsize) ){
             fposn = oposn + isrc*NBL->bsize;
             if( znzseek(fp, (long)fposn, SEEK_SET) < 0 ){
                fprintf(stderr,"** failed to locate brick %d in file '%s'\n",
                        isrc, nim->iname ? nim->iname : nim->fname);
                return -1;
             }
          }

          /* only 10,000 lines later and we're actually reading something! */
          rv = nifti_read_buffer(fp, NBL->bricks[idest], NBL->bsize, nim);
          if( rv != NBL->bsize ){
             fprintf(stderr,"** failed to read brick %d from file '%s'\n",
                     isrc, nim->iname ? nim->iname : nim->fname);
             return -1;
          }
          fposn += NBL->bsize;
       } else {
          /* we have already read this sub-brick, just copy the previous one */
          /* note that this works because they are sorted */
          memcpy(NBL->bricks[idest], NBL->bricks[sindex[c-1]], NBL->bsize);
       }

       prev = isrc;  /* in any case, note the now previous sub-brick */
   }

   return 0;
}


/*----------------------------------------------------------------------
 * nifti_alloc_NBL_mem      - allocate memory for bricks
 *
 * return 0 on success, -1 on failure
 *----------------------------------------------------------------------*/
static int nifti_alloc_NBL_mem(nifti_image * nim, int nbricks,
                               nifti_brick_list * nbl)
{
   int c;

   /* if nbricks is not specified, use the default */
   if( nbricks > 0 ) nbl->nbricks = nbricks;
   else {  /* I missed this one with the 1.17 change    02 Mar 2006 [rickr] */
      nbl->nbricks = 1;
      for( c = 4; c <= nim->ndim; c++ )
          nbl->nbricks *= nim->dim[c];
   }

   nbl->bsize   = (size_t)nim->nx * nim->ny * nim->nz * nim->nbyper;/* bytes */
   nbl->bricks  = (void **)malloc(nbl->nbricks * sizeof(void *));

   if( ! nbl->bricks ){
      fprintf(stderr,"** NANM: failed to alloc %d void ptrs\n",nbricks);
      return -1;
   }

   for( c = 0; c < nbl->nbricks; c++ ){
      nbl->bricks[c] = (void *)malloc(nbl->bsize);
      if( ! nbl->bricks[c] ){
         fprintf(stderr,"** NANM: failed to alloc %u bytes for brick %d\n",
                 (unsigned int)nbl->bsize, c);
         /* so free and clear everything before returning */
         while( c > 0 ){
            c--;
            free(nbl->bricks[c]);
         }
         free(nbl->bricks);
         nbl->bricks = NULL;
         nbl->bsize = nbl->nbricks = 0;
         return -1;
      }
   }

   if( g_opts.debug > 2 )
      fprintf(stderr,"+d NANM: alloc'd %d bricks of %u bytes for NBL\n",
              nbl->nbricks, (unsigned int)nbl->bsize);

   return 0;
}


/*----------------------------------------------------------------------
 * nifti_copynsort      - copy int list, and sort with indices
 *
 * 1. duplicate the incoming list
 * 2. create an sindex list, and init with 0..nbricks-1
 * 3. do a slow insertion sort on the small slist, along with sindex list
 * 4. check results, just to be positive
 * 
 * So slist is sorted, and sindex hold original positions.
 *
 * return 0 on success, -1 on failure
 *----------------------------------------------------------------------*/
static int nifti_copynsort(int nbricks, const int * blist, int ** slist,
                           int ** sindex)
{
   int * stmp, * itmp;   /* for ease of typing/reading */
   int   c1, c2, spos, tmp;

   *slist  = (int *)malloc(nbricks * sizeof(int));
   *sindex = (int *)malloc(nbricks * sizeof(int));

   if( !*slist || !*sindex ){
      fprintf(stderr,"** NCS: failed to alloc %d ints for sorting\n",nbricks);
      if(*slist)  free(*slist);   /* maybe one succeeded */
      if(*sindex) free(*sindex);
      return -1;
   }

   /* init the lists */
   memcpy(*slist, blist, nbricks*sizeof(int));
   for( c1 = 0; c1 < nbricks; c1++ ) (*sindex)[c1] = c1;

   /* now actually sort slist */
   stmp = *slist;
   itmp = *sindex;
   for( c1 = 0; c1 < nbricks-1; c1++ ) {
      /* find smallest value, init to current */
      spos = c1;
      for( c2 = c1+1; c2 < nbricks; c2++ )
         if( stmp[c2] < stmp[spos] ) spos = c2;
      if( spos != c1 ) /* swap: fine, don't maintain sub-order, see if I care */
      {
         tmp        = stmp[c1];      /* first swap the sorting values */
         stmp[c1]   = stmp[spos];
         stmp[spos] = tmp;

         tmp        = itmp[c1];      /* then swap the index values */
         itmp[c1]   = itmp[spos];
         itmp[spos] = tmp;
      }
   }

   if( g_opts.debug > 2 ){
      fprintf(stderr,  "+d sorted indexing list:\n");
      fprintf(stderr,  "  orig   : ");
      for( c1 = 0; c1 < nbricks; c1++ ) fprintf(stderr,"  %d",blist[c1]);
      fprintf(stderr,"\n  new    : ");
      for( c1 = 0; c1 < nbricks; c1++ ) fprintf(stderr,"  %d",stmp[c1]);
      fprintf(stderr,"\n  indices: ");
      for( c1 = 0; c1 < nbricks; c1++ ) fprintf(stderr,"  %d",itmp[c1]);
      fputc('\n', stderr);
   }

   /* check the sort (why not?  I've got time...) */
   for( c1 = 0; c1 < nbricks-1; c1++ ){
       if( (stmp[c1] > stmp[c1+1]) || (blist[itmp[c1]] != stmp[c1]) ){