fslio.c

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

    *mmy = stdmat.m[1][0] * voxx + stdmat.m[1][1] * voxy + stdmat.m[1][2] * voxz 
        + stdmat.m[1][3];
    *mmz = stdmat.m[2][0] * voxx + stdmat.m[2][1] * voxy + stdmat.m[2][2] * voxz 
        + stdmat.m[2][3];
}


void FslGetVoxCoord(mat44 stdmat, float mmx, float mmy, float mmz, 
                   float *voxx, float *voxy, float *voxz) 
{
  mat44 mm2vox;

  mm2vox = nifti_mat44_inverse(stdmat);
    *voxx = mm2vox.m[0][0] * mmx + mm2vox.m[0][1] * mmy + mm2vox.m[0][2] * mmz 
        + mm2vox.m[0][3];
    *voxy = mm2vox.m[1][0] * mmx + mm2vox.m[1][1] * mmy + mm2vox.m[1][2] * mmz 
        + mm2vox.m[1][3];
    *voxz = mm2vox.m[2][0] * mmx + mm2vox.m[2][1] * mmy + mm2vox.m[2][2] * mmz 
        + mm2vox.m[2][3];
}


void FslSetStdXform(FSLIO *fslio, short sform_code, mat44 stdmat)
{
    /* NB: stdmat must point to a 4x4 array */
  if (fslio==NULL)  FSLIOERR("FslSetStdXform: Null pointer passed for FSLIO");
  if (fslio->niftiptr!=NULL) {
      fslio->niftiptr->sform_code = sform_code;
      fslio->niftiptr->sto_xyz.m[0][0] = stdmat.m[0][0];
      fslio->niftiptr->sto_xyz.m[0][1] = stdmat.m[0][1];
      fslio->niftiptr->sto_xyz.m[0][2] = stdmat.m[0][2];
      fslio->niftiptr->sto_xyz.m[0][3] = stdmat.m[0][3];
      fslio->niftiptr->sto_xyz.m[1][0] = stdmat.m[1][0];
      fslio->niftiptr->sto_xyz.m[1][1] = stdmat.m[1][1];
      fslio->niftiptr->sto_xyz.m[1][2] = stdmat.m[1][2];
      fslio->niftiptr->sto_xyz.m[1][3] = stdmat.m[1][3];
      fslio->niftiptr->sto_xyz.m[2][0] = stdmat.m[2][0];
      fslio->niftiptr->sto_xyz.m[2][1] = stdmat.m[2][1];
      fslio->niftiptr->sto_xyz.m[2][2] = stdmat.m[2][2];
      fslio->niftiptr->sto_xyz.m[2][3] = stdmat.m[2][3];
      fslio->niftiptr->sto_xyz.m[3][0] = 0;
      fslio->niftiptr->sto_xyz.m[3][1] = 0;
      fslio->niftiptr->sto_xyz.m[3][2] = 0;
      fslio->niftiptr->sto_xyz.m[3][3] = 1;
      fslio->niftiptr->sto_ijk = nifti_mat44_inverse(fslio->niftiptr->sto_xyz);
  }
  if (fslio->mincptr!=NULL) {
    fprintf(stderr,"Warning:: Minc is not yet supported\n");
  }
}


short FslGetStdXform(FSLIO *fslio, mat44 *stdmat)
{
    /* returns sform code  (NB: stdmat must point to a 4x4 array) */
    float dx,dy,dz,tr;
    if (fslio==NULL)  FSLIOERR("FslGetStdXform: Null pointer passed for FSLIO");
    if (fslio->niftiptr!=NULL) {
        stdmat->m[0][0] = fslio->niftiptr->sto_xyz.m[0][0];
        stdmat->m[0][1] = fslio->niftiptr->sto_xyz.m[0][1];
        stdmat->m[0][2] = fslio->niftiptr->sto_xyz.m[0][2];
        stdmat->m[0][3] = fslio->niftiptr->sto_xyz.m[0][3];
        stdmat->m[1][0] = fslio->niftiptr->sto_xyz.m[1][0];
        stdmat->m[1][1] = fslio->niftiptr->sto_xyz.m[1][1];
        stdmat->m[1][2] = fslio->niftiptr->sto_xyz.m[1][2];
        stdmat->m[1][3] = fslio->niftiptr->sto_xyz.m[1][3];
        stdmat->m[2][0] = fslio->niftiptr->sto_xyz.m[2][0];
        stdmat->m[2][1] = fslio->niftiptr->sto_xyz.m[2][1];
        stdmat->m[2][2] = fslio->niftiptr->sto_xyz.m[2][2];
        stdmat->m[2][3] = fslio->niftiptr->sto_xyz.m[2][3];
        stdmat->m[3][0] = 0.0;
        stdmat->m[3][1] = 0.0;
        stdmat->m[3][2] = 0.0;
        stdmat->m[3][3] = 1.0;
        
        /* the code below gives a default but it really should never be used */
        if (fslio->niftiptr->sform_code == NIFTI_XFORM_UNKNOWN) {
            FslGetVoxDim(fslio,&dx,&dy,&dz,&tr);
            stdmat->m[0][0] = -dx;  /* default Radiological convention */
            stdmat->m[0][1] = 0;
            stdmat->m[0][2] = 0;
            stdmat->m[0][3] = 0;
            stdmat->m[1][0] = 0;
            stdmat->m[1][1] = dy;
            stdmat->m[1][2] = 0;
            stdmat->m[1][3] = 0;
            stdmat->m[2][0] = 0;
            stdmat->m[2][1] = 0;
            stdmat->m[2][2] = dz;
            stdmat->m[2][3] = 0;
            stdmat->m[3][0] = 0.0;
            stdmat->m[3][1] = 0.0;
            stdmat->m[3][2] = 0.0;
            stdmat->m[3][3] = 1.0;
        }
        return fslio->niftiptr->sform_code;
    }
    if (fslio->mincptr!=NULL) {
        fprintf(stderr,"Warning:: Minc is not yet supported\n");
    }
    return NIFTI_XFORM_UNKNOWN;
}


void FslSetRigidXform(FSLIO *fslio, short qform_code, mat44 rigidmat)
{
    /* NB: rigidmat must point to an allocated mat44 */
  float dx, dy, dz;
  if (fslio==NULL)  FSLIOERR("FslSetRigidXform: Null pointer passed for FSLIO");
  if (fslio->niftiptr!=NULL) {
      fslio->niftiptr->qform_code = qform_code;
      fslio->niftiptr->qto_xyz.m[0][0] = rigidmat.m[0][0];
      fslio->niftiptr->qto_xyz.m[0][1] = rigidmat.m[0][1];
      fslio->niftiptr->qto_xyz.m[0][2] = rigidmat.m[0][2];
      fslio->niftiptr->qto_xyz.m[0][3] = rigidmat.m[0][3];
      fslio->niftiptr->qto_xyz.m[1][0] = rigidmat.m[1][0];
      fslio->niftiptr->qto_xyz.m[1][1] = rigidmat.m[1][1];
      fslio->niftiptr->qto_xyz.m[1][2] = rigidmat.m[1][2];
      fslio->niftiptr->qto_xyz.m[1][3] = rigidmat.m[1][3];
      fslio->niftiptr->qto_xyz.m[2][0] = rigidmat.m[2][0];
      fslio->niftiptr->qto_xyz.m[2][1] = rigidmat.m[2][1];
      fslio->niftiptr->qto_xyz.m[2][2] = rigidmat.m[2][2];
      fslio->niftiptr->qto_xyz.m[2][3] = rigidmat.m[2][3];
      fslio->niftiptr->qto_xyz.m[3][0] = 0;
      fslio->niftiptr->qto_xyz.m[3][1] = 0;
      fslio->niftiptr->qto_xyz.m[3][2] = 0;
      fslio->niftiptr->qto_xyz.m[3][3] = 1;
      nifti_mat44_to_quatern(
            fslio->niftiptr->qto_xyz,&(fslio->niftiptr->quatern_b),
            &(fslio->niftiptr->quatern_c),&(fslio->niftiptr->quatern_d),
            &(fslio->niftiptr->qoffset_x),&(fslio->niftiptr->qoffset_y),
            &(fslio->niftiptr->qoffset_z),&dx,&dy,&dz,&(fslio->niftiptr->qfac));
      fslio->niftiptr->qto_ijk = nifti_mat44_inverse(fslio->niftiptr->qto_xyz);

  }
  if (fslio->mincptr!=NULL) {
    fprintf(stderr,"Warning:: Minc is not yet supported\n");
  }
}


short FslGetRigidXform(FSLIO *fslio, mat44 *rigidmat)
{
    /* returns qform code  (NB: rigidmat must point to an allocated mat44) */
    float dx,dy,dz,tr;
    if (fslio==NULL)  FSLIOERR("FslGetRigidXform: Null pointer passed for FSLIO");
    if (fslio->niftiptr!=NULL) {
        rigidmat->m[0][0] = fslio->niftiptr->qto_xyz.m[0][0];
        rigidmat->m[0][1] = fslio->niftiptr->qto_xyz.m[0][1];
        rigidmat->m[0][2] = fslio->niftiptr->qto_xyz.m[0][2];
        rigidmat->m[0][3] = fslio->niftiptr->qto_xyz.m[0][3];
        rigidmat->m[1][0] = fslio->niftiptr->qto_xyz.m[1][0];
        rigidmat->m[1][1] = fslio->niftiptr->qto_xyz.m[1][1];
        rigidmat->m[1][2] = fslio->niftiptr->qto_xyz.m[1][2];
        rigidmat->m[1][3] = fslio->niftiptr->qto_xyz.m[1][3];
        rigidmat->m[2][0] = fslio->niftiptr->qto_xyz.m[2][0];
        rigidmat->m[2][1] = fslio->niftiptr->qto_xyz.m[2][1];
        rigidmat->m[2][2] = fslio->niftiptr->qto_xyz.m[2][2];
        rigidmat->m[2][3] = fslio->niftiptr->qto_xyz.m[2][3];
        rigidmat->m[3][0] = 0.0;
        rigidmat->m[3][1] = 0.0;
        rigidmat->m[3][2] = 0.0;
        rigidmat->m[3][3] = 1.0;
        
        /* the code gives a default but it should never really be used */
        if (fslio->niftiptr->sform_code == NIFTI_XFORM_UNKNOWN) {
          FslGetVoxDim(fslio,&dx,&dy,&dz,&tr);
          rigidmat->m[0][0] = dx;
          rigidmat->m[0][1] = 0;
          rigidmat->m[0][2] = 0;
          rigidmat->m[0][3] = 0;
          rigidmat->m[1][0] = 0;
          rigidmat->m[1][1] = dy;
          rigidmat->m[1][2] = 0;
          rigidmat->m[1][3] = 0;
          rigidmat->m[2][0] = 0;
          rigidmat->m[2][1] = 0;
          rigidmat->m[2][2] = dz;
          rigidmat->m[3][0] = 0.0;
          rigidmat->m[3][1] = 0.0;
          rigidmat->m[3][2] = 0.0;
          rigidmat->m[3][3] = 1.0;
        }
        return fslio->niftiptr->qform_code;
    }
    if (fslio->mincptr!=NULL) {
        fprintf(stderr,"Warning:: Minc is not yet supported\n");
    }
    return NIFTI_XFORM_UNKNOWN;
}


void FslSetIntent(FSLIO *fslio, short intent_code, float p1, float p2, float p3)
{
  if (fslio==NULL)  FSLIOERR("FslSetIntent: Null pointer passed for FSLIO");
  if (fslio->niftiptr!=NULL) {
      fslio->niftiptr->intent_code = intent_code;
      fslio->niftiptr->intent_p1 = p1;
      fslio->niftiptr->intent_p2 = p2;
      fslio->niftiptr->intent_p3 = p3;
  }
  if (fslio->mincptr!=NULL) {
    fprintf(stderr,"Warning:: Minc is not yet supported\n");
  }
}


short FslGetIntent(FSLIO *fslio, short *intent_code, float *p1, float *p2,
                   float *p3)
{
  /* also returns intent code */
  if (fslio==NULL)  FSLIOERR("FslGetIntent: Null pointer passed for FSLIO");
  if (fslio->niftiptr!=NULL) {
      *intent_code = fslio->niftiptr->intent_code;
      *p1 = fslio->niftiptr->intent_p1;
      *p2 = fslio->niftiptr->intent_p2;
      *p3 = fslio->niftiptr->intent_p3;
      return fslio->niftiptr->intent_code;
  }
  if (fslio->mincptr!=NULL) {
    fprintf(stderr,"Warning:: Minc is not yet supported\n");
  }
  return NIFTI_INTENT_NONE;
}




void FslSetIntensityScaling(FSLIO *fslio, float slope, float intercept)
{
  if (fslio==NULL)  FSLIOERR("FslSetIntensityScaling: Null pointer passed for FSLIO");
  if (fslio->niftiptr!=NULL) {
      fslio->niftiptr->scl_slope = slope;
      fslio->niftiptr->scl_inter = intercept;
  }
  if (fslio->mincptr!=NULL) {
    fprintf(stderr,"Warning:: Minc is not yet supported\n");
  }
}


int FslGetIntensityScaling(FSLIO *fslio, float *slope, float *intercept)
{
  /* returns 1 if scaling required or 0 otherwise */
  if (fslio==NULL)  FSLIOERR("FslGetIntensityScaling: Null pointer passed for FSLIO");
  if (fslio->niftiptr!=NULL) {
    *slope = fslio->niftiptr->scl_slope;
    *intercept = fslio->niftiptr->scl_inter;
    if (fabs(*slope)<1e-30) {
      *slope = 1.0;
      *intercept = 0.0;
      return 0;
    }
    if ( (fabs(*slope - 1.0)>1e-30) || (fabs(*intercept)>1e-30) ) {
      return 1;
    } else {
      return 0;
    }
  }
  if (fslio->mincptr!=NULL) {
    fprintf(stderr,"Warning:: Minc is not yet supported\n");
  }
  return 0;
 
}


mat33 mat44_to_mat33(mat44 x)
{
  mat33 y;
  int i,j;
  for (i=0; i<3; i++) {
    for (j=0; j<3; j++) {
      y.m[i][j] = x.m[i][j];
    }
  }
  return y;
}


int FslGetLeftRightOrder(FSLIO *fslio)
{
  /* Determines if the image is stored in neurological or radiological convention */
  int order=FSL_RADIOLOGICAL, sform_code, qform_code;
  float det=-1.0;
  mat44 sform44, qform44;
  mat33 sform33, qform33;
  if (fslio==NULL)  FSLIOERR("FslGetLeftRightOrder: Null pointer passed for FSLIO");
  if (fslio->niftiptr!=NULL) {
    sform_code = FslGetStdXform(fslio,&sform44);
    qform_code = FslGetRigidXform(fslio,&qform44);
    if (sform_code!=NIFTI_XFORM_UNKNOWN) { 
      sform33 = mat44_to_mat33(sform44);
      det = nifti_mat33_determ(sform33);
    } else if (qform_code!=NIFTI_XFORM_UNKNOWN) { 
      qform33 = mat44_to_mat33(qform44);
      det = nifti_mat33_determ(qform33); 
    }

    if (det<0.0) order=FSL_RADIOLOGICAL;
    else order=FSL_NEUROLOGICAL;
  }
  if (fslio->mincptr!=NULL) {
    fprintf(stderr,"Warning:: Minc is not yet supported\n");
  }
  return order;
}



void FslSetAnalyzeSform(FSLIO *fslio, const short *orig,
                        float dx, float dy, float dz)
{
  /* Creates an sform matrix for an Analyze file */
  /* THIS ALWAYS CREATES A RADIOLOGICAL ORDERED SFORM */
  /* NB: the origin passed in here is in Analyze convention - starting at 1, not 0 */
  float x, y, z;
  if (fslio==NULL)  FSLIOERR("FslSetAnalyzeSform: Null pointer passed for FSLIO");
  if (fslio->niftiptr!=NULL) {
    if (FslBaseFileType(FslGetFileType(fslio))==FSL_TYPE_ANALYZE) {
      /* default case */
      fslio->niftiptr->sform_code = NIFTI_XFORM_UNKNOWN;
    }
    /* ignore all zero origins - really all serious coord stuff should
       be done via the FslSetStdCoord call */
    if ((orig[0]!=0) || (orig[1]!=0) || (orig[2]!=0))
      {
        short origx=0, origy=0, origz=0;
        if ((orig[0]!=0) || (orig[1]!=0) || (orig[2]!=0)) {
          /* convert to nifti conventions (start at 0 not 1) */
          origx = orig[0] - 1;
          origy = orig[1] - 1;
          origz = orig[2] - 1;
        }
        if ( dx * dy * dz > 0 ) {
          /* change neurological convention to radiological if necessary */
          dx = -dx;
        }
        if ( (FslBaseFileType(FslGetFileType(fslio))==FSL_TYPE_ANALYZE) 
             || (fslio->niftiptr->sform_code == NIFTI_XFORM_UNKNOWN) ) {
          /* make a default transform with the requested origin at xyz=000 */ 
          fslio->niftiptr->sform_code = NIFTI_XFORM_ALIGNED_ANAT;
          fslio->niftiptr->sto_xyz.m[0][0] = dx;
          fslio->niftiptr->sto_xyz.m[0][1] = 0;
          fslio->niftiptr->sto_xyz.m[0][2] = 0;
          fslio->niftiptr->sto_xyz.m[0][3] = -(origx)*(dx);
          fslio->niftiptr->sto_xyz.m[1][0] = 0;
          fslio->niftiptr->sto_xyz.m[1][1] = dy;
          fslio->niftiptr->sto_xyz.m[1][2] = 0;
          fslio->niftiptr->sto_xyz.m[1][3] = -(origy)*(dy);
          fslio->niftiptr->sto_xyz.m[2][0] = 0;
          fslio->niftiptr->sto_xyz.m[2][1] = 0;
          fslio->niftiptr->sto_xyz.m[2][2] = dz;
          fslio->niftiptr->sto_xyz.m[2][3] = -(origz)*(dz);
          fslio->niftiptr->sto_xyz.m[3][0] = 0;
          fslio->niftiptr->sto_xyz.m[3][1] = 0;
          fslio->niftiptr->sto_xyz.m[3][2] = 0;
          fslio->niftiptr->sto_xyz.m[3][3] = 1;
          fslio->niftiptr->sto_ijk =
                 nifti_mat44_inverse(fslio->niftiptr->sto_xyz);
        } else {
          /* update the existing origin */
          /* find out what the existing xyz of the requested origin is */
          x = fslio->niftiptr->sto_xyz.m[0][0] * origx
            + fslio->niftiptr->sto_xyz.m[0][1] * origy