dcmimage.c

Go to the documentation of this file.

 
/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <string.h>
/*****************************************************************************/
#include "tpcextensions.h"
/*****************************************************************************/
#include "tpcdcm.h"
/*****************************************************************************/
 
/*****************************************************************************/
int dcmImgIsotope(
  DCMFILE *d,
  isotope *isot,
  decaycorrection *dc,
  const int verbose
) {
  if(verbose>0) {printf("%s(dcmfile, *isotope, *decaycorrection)\n", __func__); fflush(stdout);}
  if(isot!=NULL) *isot=ISOTOPE_UNKNOWN;
  if(dc!=NULL) *dc=DECAY_UNKNOWN;
  if(d==NULL || d->item==NULL) return(1);
 
 
  DCMTAG tag; DCMITEM *iptr, *jptr;
 
  /* Get the status of decay correction */
  if(dc!=NULL) {
    if(verbose>1) {printf("trying to find 'Decay Correction' tag\n"); fflush(stdout);}
    dcmTagSet(&tag, 0x0054, 0x1102); iptr=dcmFindTag(d->item, 0, &tag, 0);
    if(iptr!=NULL) {
      char *buf=dcmValueString(iptr); 
      if(verbose>12) printf("Decay Correction := '%s'\n", buf);
      if(strcasestr(buf, "NONE")!=NULL) *dc=DECAY_NOTCORRECTED;
      else if(strcasestr(buf, "START")!=NULL) *dc=DECAY_CORRECTED_START;
      else if(strcasestr(buf, "ADMIN")!=NULL) *dc=DECAY_CORRECTED_ADMIN;
      free(buf);
    } else {
      if(verbose>1) {printf("  not found\ntrying to find 'Decay Corrected' tag\n"); fflush(stdout);}
      dcmTagSet(&tag, 0x0018, 0x9758); iptr=dcmFindTag(d->item, 0, &tag, 0);
      if(iptr!=NULL) {
        char *buf=dcmValueString(iptr); 
        if(verbose>12) printf("Decay Correction := '%s'\n", buf);
        if(strcasestr(buf, "YES")!=NULL) *dc=DECAY_CORRECTED;
        else if(strcasestr(buf, "NO")!=NULL) *dc=DECAY_NOTCORRECTED;
        free(buf);
      } else if(verbose>1) {printf("  not found.\n"); fflush(stdout);}
    }
    if(verbose>0 && *dc==DECAY_UNKNOWN) {
      fprintf(stderr, "Warning: status of decay correction not found.\n");
      fflush(stdout);
    }
    if(verbose>11) printf("decayCorrection := %s\n", decayDescr(*dc));
  }
 
  if(isot==NULL) return(0);
 
  /* Get the radioisotope */
  if(verbose>1) {printf("trying to find Radiopharmaceutical Information Sequence\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x0054, 0x0016); iptr=dcmFindTag(d->item, 0, &tag, 0);
  if(iptr!=NULL) {
    if(verbose>1) {printf("  found.\n"); fflush(stdout);}
    iptr=iptr->child_item;
  } else {
    if(verbose>1) {printf("  not found.\n"); fflush(stdout);}
    iptr=d->item;
  }
  if(verbose>1) {printf("trying to find Radionuclide Code Sequence\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x0054, 0x0300); jptr=dcmFindTag(iptr, 0, &tag, 0);
  if(jptr!=NULL) {
    if(verbose>1) {printf("  found.\n"); fflush(stdout);}
    if(verbose>1) {printf("trying to find Code Meaning\n"); fflush(stdout);}
    dcmTagSet(&tag, 0x0008, 0x0104); jptr=dcmFindDownTag(jptr, 0, &tag, 0);
    if(jptr!=NULL) {
      char *buf=dcmValueString(jptr); 
      if(verbose>12) {printf("Code Meaning := '%s'\n", buf); fflush(stdout);}
      *isot=isotopeIdentify(buf);
      free(buf);
      if(verbose>11) {printf("isotope := %s\n", isotopeName(*isot)); fflush(stdout);}
    } else if(verbose>1) {printf("  not found.\n"); fflush(stdout);}
  } else {
    if(verbose>1) {printf("  not found.\n"); fflush(stdout);}
  }
  if(*isot==ISOTOPE_UNKNOWN) {
    if(verbose>1) {printf("trying to find Radionuclide Half Life\n"); fflush(stdout);}
    dcmTagSet(&tag, 0x0018, 0x1075); jptr=dcmFindTag(iptr, 0, &tag, 0);
    if(jptr!=NULL) {
      char *buf=dcmValueString(jptr); strCleanSpaces(buf);
      if(verbose>12) {printf("Radionuclide Half Life := '%s'\n", buf); fflush(stdout);}
      double hl=atofVerified(buf);
      *isot=isotopeIdentifyHalflife(hl/60.0);
      free(buf);
      if(verbose>11) {printf("isotope := %s\n", isotopeName(*isot)); fflush(stdout);}
    } else if(verbose>1) {printf("  not found.\n"); fflush(stdout);}
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int dcmImgPos(
  DCMFILE *d,
  double *imgpos,
  const int verbose
) {
  if(verbose>0) {printf("%s(dcmfile, imgpos)\n", __func__); fflush(stdout);}
  if(d==NULL || d->item==NULL || imgpos==NULL) return(1);
  for(int i=0; i<3; i++) imgpos[i]=nan("");
 
  DCMTAG tag; DCMITEM *iptr, *jptr;
 
  if(verbose>1) {printf("trying to find Per Frame Functional Groups Sequence\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x5200, 0x9230); // Per Frame Functional Groups Sequence
  iptr=dcmFindTag(d->item, 0, &tag, 0);
  if(iptr!=NULL) {
    if(verbose>1) {printf("  found.\n"); fflush(stdout);}
 
    /* search for Frame Content Sequences until we find one with position number 1 */
    int posnr=0; dcmTagSet(&tag, 0x0020, 0x9111);
    jptr=dcmFindDownTag(iptr->child_item, 0, &tag, 0);
    while(jptr!=NULL) {
      /* Get In Stack Position Number */
      DCMITEM *kptr; DCMTAG spostag; dcmTagSet(&spostag, 0x0020, 0x9057);
      kptr=dcmFindDownTag(jptr->child_item, 0, &spostag, 0);
      if(kptr!=NULL) {
        posnr=dcmitemGetInt(kptr); if(verbose>3) printf("  posnr=%d\n", posnr);
        if(posnr==1) break; // found stack position number 1
      }
      jptr=dcmFindDownTag(jptr->next_item, 0, &tag, 0);
    }
    /* If stack position 1 was found, then get plane position sequence */
    if(posnr==1 && jptr!=NULL) {
      dcmTagSet(&tag, 0x0020, 0x9113);
      jptr=dcmFindDownTag(jptr->next_item, 0, &tag, 0);
      if(jptr!=NULL) {
        /* Get Image Position (Patient) under it */
        DCMITEM *kptr; DCMTAG ipostag; dcmTagSet(&ipostag, 0x0020, 0x0032);
        kptr=dcmFindDownTag(jptr->child_item, 0, &ipostag, 0);
        if(kptr!=NULL) {
          char *buf=dcmValueString(kptr);
          if(verbose>4) printf("  image_position := '%s'\n", buf);
          int n=sscanf(buf, "%lf\\%lf\\%lf", &imgpos[0], &imgpos[1], &imgpos[2]);
          free(buf);
          if(n!=3) {
            if(verbose>0) fprintf(stderr, "Error: invalid Image Position (Patient)\n");
            return(131);
          }
        } else {
          if(verbose>0) fprintf(stderr, "Error: no Image Position (Patient)\n");
          return(121);
        }
      } else {
        if(verbose>0) fprintf(stderr, "Error: no Plane Position Sequence\n");
        return(111);
      }
    } else {
      if(verbose>0) {fprintf(stderr, "Error: stack position 1 not found.\n"); fflush(stderr);}
      return(101);
    }
    return(0);
  }
  if(verbose>1) {printf("  not found.\n"); fflush(stdout);}
 
  if(verbose>1) {printf("find the smallest Instance (image) Number\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x0020, 0x0013);
  iptr=dcmFindTag(d->item, 0, &tag, 0);
  if(iptr==NULL) {
    if(verbose>0) {fprintf(stderr, "Error: Instance (image) Number not found.\n"); fflush(stderr);}
    return(201);
  }
  int mininr=dcmitemGetInt(iptr); if(verbose>5) printf("  inr=%d\n", mininr);
  jptr=iptr;
  while(jptr!=NULL) {
    jptr=dcmFindDownTag(jptr->next_item, 0, &tag, 0);
    if(jptr==NULL) break;
    int inr=dcmitemGetInt(jptr); if(verbose>5) printf("  inr=%d\n", inr);
    if(inr<mininr) {mininr=inr; iptr=jptr;}
  }
  if(verbose>2) {printf("  smallest Instance (image) Number: %d\n", mininr); fflush(stdout);}
 
  /* Get Image Position (Patient) after the smallest Instance Number */
  dcmTagSet(&tag, 0x0020, 0x0032);
  jptr=dcmFindDownTag(iptr, 0, &tag, 0);
  if(jptr==NULL) {
    if(verbose>0) fprintf(stderr, "Error: no Image Position (Patient)\n");
    return(211);
  } else {
    char *buf=dcmValueString(jptr);
    if(verbose>2) printf("  image_position := '%s'\n", buf);
    int n=sscanf(buf, "%lf\\%lf\\%lf", &imgpos[0], &imgpos[1], &imgpos[2]);
    free(buf);
    if(n!=3) {
      if(verbose>0) fprintf(stderr, "Error: invalid Image Position (Patient)\n");
      return(212);
    }
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int dcmImgDim(
  DCMFILE *d,
  unsigned short int *imgdim,
  const int verbose
) {
  if(verbose>0) {printf("%s(dcmfile, imgdim)\n", __func__); fflush(stdout);}
  if(d==NULL || d->item==NULL || imgdim==NULL) return(1);
  for(int i=0; i<4; i++) imgdim[i]=0;
 
  DCMTAG tag; DCMITEM *iptr;
 
  /* Get image row number (dimy) */
  dcmTagSet(&tag, 0x0028, 0x0010);
  iptr=dcmFindTag(d->item, 0, &tag, 0);
  if(iptr==NULL) {
    if(verbose>0) {fprintf(stderr, "Error: cannot find matrix row number.\n"); fflush(stderr);}
    return(11);
  }
  imgdim[1]=(unsigned short int)dcmitemGetInt(iptr);
 
  /* Get image column number (dimx) */
  dcmTagSet(&tag, 0x0028, 0x0011);
  iptr=dcmFindTag(d->item, 0, &tag, 0);
  if(iptr==NULL) {
    if(verbose>0) {fprintf(stderr, "Error: cannot find matrix column number.\n"); fflush(stderr);}
    return(21);
  }
  imgdim[0]=(unsigned short int)dcmitemGetInt(iptr);
 
  /* Get image plane/slice number (dimz) */
  int a, b;
  if(verbose>1) {printf("trying to read plane number as 'Number of Slices'\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x0054, 0x0081);
  iptr=dcmFindTag(d->item, 0, &tag, 0);
  if(iptr!=NULL) {
    imgdim[2]=(unsigned short int)dcmitemGetInt(iptr);
    if(verbose>2) printf("planeNr := %u\n", imgdim[2]); 
  } else {
    if(verbose>1) {printf("trying to read plane number from 'In Stack Position Number'\n"); fflush(stdout);}
    dcmTagSet(&tag, 0x0020, 0x9057);
    if(dcmTagIntRange(d->item, &tag, &a, &b, verbose-100)==0) {
      imgdim[2]=(unsigned short int)b; 
      if(verbose>2) printf("planeNr := %u\n", imgdim[2]); 
    } else {
      if(verbose>1) {printf("trying to read plane number from 'Dimension Index Values'\n"); fflush(stdout);}
      dcmTagSet(&tag, 0x0020, 0x9157);
      if(dcmTagIntRange(d->item, &tag, &a, &b, verbose-1)==0) {
        imgdim[2]=(unsigned short int)b; 
        if(verbose>2) printf("planeNr := %u\n", imgdim[2]); 
      } else {
        if(verbose>0) {fprintf(stderr, "Warning: cannot find plane number.\n"); fflush(stderr);}
        imgdim[2]=1;
      }
    }
  }
 
 
  /* Get frame number */
  if(verbose>1) {printf("trying to get frame number from 'Temporal position index'\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x0020, 0x9128);
  if(dcmTagIntRange(d->item, &tag, &a, &b, verbose-100)==0) {
    imgdim[3]=(unsigned short int)b; 
    if(verbose>2) printf("frameNr := %u\n", imgdim[3]); 
  } else {
    if(verbose>1) {printf("trying to get frame number from 'Number of Time Slices'\n"); fflush(stdout);}
    dcmTagSet(&tag, 0x0054, 0x0101);
    iptr=dcmFindTag(d->item, 0, &tag, 0);
    if(iptr!=NULL) {
      imgdim[3]=(unsigned short int)dcmitemGetInt(iptr);
      if(verbose>2) printf("frameNr := %u\n", imgdim[3]);
    } else {
      if(verbose>0) {fprintf(stderr, "Warning: cannot find frame number.\n"); fflush(stderr);}
      imgdim[3]=1;
    }
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int dcmImgPxlsize(
  DCMFILE *d,
  double *pxlsize,
  const int verbose
) {
  if(verbose>0) {printf("%s(dcmfile, pxlsize)\n", __func__); fflush(stdout);}
  if(d==NULL || d->item==NULL || pxlsize==NULL) return(1);
  for(int i=0; i<3; i++) pxlsize[i]=nan("");
 
  DCMTAG tag; DCMITEM *iptr, *jptr, *kptr;
 
  if(verbose>1) {printf("trying to find Pixel Measures Sequence\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x0028, 0x9110);
  iptr=dcmFindTag(d->item, 0, &tag, 0);
  if(iptr==NULL) {
    kptr=d->item; // search from root
    if(verbose>3) {fprintf(stdout, "Note: missing Pixel Measures Sequence.\n"); fflush(stdout);}
  } else {
    kptr=iptr->child_item; // search under this sequence
  }
  dcmTagSet(&tag, 0x0018, 0x0050);
  jptr=dcmFindDownTag(kptr, 0, &tag, 0);
  if(jptr==NULL) {
    if(verbose>0) {fprintf(stderr, "Error: missing Slice Thickness.\n"); fflush(stderr);}
    return(11);
  }
  pxlsize[2]=dcmitemGetReal(jptr);
 
  if(verbose>1) {printf("trying to find Spacing Between Slices\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x0018, 0x0088);
  jptr=dcmFindDownTag(kptr, 0, &tag, 0);
  if(jptr==NULL) {
    if(verbose>1) {fprintf(stderr, "Warning: missing Spacing Between Slices.\n"); fflush(stderr);}
  } else {
    double v=dcmitemGetReal(jptr);
    if(v>pxlsize[2]) pxlsize[2]=v;
  }
 
  if(verbose>1) {printf("trying to find Pixel Spacing/Size\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x0028, 0x0030);
  jptr=dcmFindDownTag(kptr, 0, &tag, 0);
  if(jptr==NULL) {
    if(verbose>0) {fprintf(stderr, "Error: missing Pixel Spacing/Size.\n"); fflush(stderr);}
    return(21);
  }
  char *buf=dcmValueString(jptr);
  int n=sscanf(buf, "%lf\\%lf", &pxlsize[1], &pxlsize[0]);
  free(buf);
  if(n!=2) {
    if(verbose>0) {fprintf(stderr, "Error: invalid Pixel Spacing/Size\n"); fflush(stderr);}
    return(22);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int dcmImgOrient(
  DCMFILE *d,
  double *iop,
  const int verbose
) {
  if(verbose>0) {printf("%s(dcmfile, iop)\n", __func__); fflush(stdout);}
  if(d==NULL || d->item==NULL || iop==NULL) return(1);
  for(int i=0; i<6; i++) iop[i]=nan("");
 
  DCMTAG tag; DCMITEM *iptr;
 
  if(verbose>1) {printf("reading Image Orientation (Patient)\n"); fflush(stdout);}
  dcmTagSet(&tag, 0x0020, 0x0037);
  iptr=dcmFindTag(d->item, 0, &tag, 0);
  if(iptr==NULL) {
    if(verbose>0) {fprintf(stderr, "Error: no Image Orientation (Patient)\n"); fflush(stderr);}
    return(11);
  }
  char *buf=dcmValueString(iptr);
  if(verbose>2) printf("  image_orientation := '%s'\n", buf);
  int n=sscanf(buf, "%lf\\%lf\\%lf\\%lf\\%lf\\%lf", &iop[0], &iop[1], &iop[2], 
                                                    &iop[3], &iop[4], &iop[5]);
  free(buf);
  if(n!=6) {
    if(verbose>0) {fprintf(stderr, "Error: invalid Image Orientation (Patient)\n"); fflush(stderr);}
    return(12);
  }
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int dcmImgXform(
  double *iop,
  double *xyzMM,
  double *imgPos,
  double *xform,
  const int verbose
) {
  if(verbose>0) {printf("%s(iop[], xyzMM[], imgPos[], xform)\n", __func__); fflush(stdout);}
  if(iop==NULL || xyzMM==NULL || imgPos==NULL || xform==NULL) return(1);
  for(int i=0; i<16; i++) xform[i]=nan("");
 
 
  double readV[3], phaseV[3];
  readV[0]=iop[0]; readV[1]=iop[1]; readV[2]=iop[2];
  phaseV[0]=iop[3]; phaseV[1]=iop[4]; phaseV[2]=iop[5];
  if(verbose>3) {
    printf("readV := %g", readV[0]); for(int i=1; i<3; i++) printf(", %g", readV[i]);
    printf("\n");
    printf("phaseV := %g", phaseV[0]); for(int i=1; i<3; i++) printf(", %g", phaseV[i]);
    printf("\n");
  }
 
  /* Cross-product of vectors */
  double sliceV[3];
  sliceV[0]= readV[1]*phaseV[2] - phaseV[1]*readV[2];
  sliceV[1]=-readV[0]*phaseV[2] + phaseV[0]*readV[2];
  sliceV[2]= readV[0]*phaseV[1] - phaseV[0]*readV[1];
  if(verbose>3) {
    printf("sliceV := %g", sliceV[0]); for(int i=1; i<3; i++) printf(", %g", sliceV[i]); 
    printf("\n");
  }
 
  /* Copy vectors into a 3x3 matrix */
  double rM[3][3];
  for(int i=0; i<3; i++) rM[0][i]=readV[i];
  for(int i=0; i<3; i++) rM[1][i]=phaseV[i];
  for(int i=0; i<3; i++) rM[2][i]=sliceV[i];
  /* Transpose the matrix */
  {
    double tM[3][3];
    for(int i=0; i<3; i++) for(int j=0; j<3; j++) tM[i][j]=rM[i][j];
    for(int i=0; i<3; i++) for(int j=0; j<3; j++) rM[i][j]=tM[j][i];
  }
  if(verbose>3) {
    printf("rM :=\t%g", rM[0][0]); 
    for(int i=1; i<3; i++) printf("\t%g", rM[0][i]);
    printf("\n");
    for(int i=0; i<3; i++) printf("\t%g", rM[1][i]);
    printf("\n");
    for(int i=0; i<3; i++) printf("\t%g", rM[2][i]);
    printf("\n");
  }
  /* Max (left-over) indices */
  unsigned short int xyz[3];
  xyz[0]=0; 
  if(fabs(rM[1][0])>fabs(rM[0][0]) && fabs(rM[1][0])>fabs(rM[2][0])) xyz[0]=1;
  else if(fabs(rM[2][0])>fabs(rM[0][0]) && fabs(rM[1][0])>fabs(rM[1][0])) xyz[0]=2;
  if(xyz[0]==0)      {if(fabs(rM[2][1])>fabs(rM[1][1])) xyz[1]=2; else xyz[1]=1;}
  else if(xyz[0]==1) {if(fabs(rM[2][1])>fabs(rM[0][1])) xyz[1]=2; else xyz[1]=0;}
  else               {if(fabs(rM[1][1])>fabs(rM[0][1])) xyz[1]=1; else xyz[1]=0;}
  xyz[2]=3-xyz[0]-xyz[1];
  if(verbose>3) {
    printf("xyz := %u", xyz[0]); for(int i=1; i<3; i++) printf(", %u", xyz[i]); 
    printf("\n(0, 1, 2 for Sag/Cor/Tra slice)\n");
  }
  double cosSL=rM[xyz[2]][2];
  if(verbose>3) {printf("cosSL := %g\n", cosSL); fflush(stdout);}
 
  /* Multiply matrix with voxel sizes */
  if(verbose>3) {
    printf("xyzMM := %g", xyzMM[0]); for(int i=1; i<3; i++) printf(", %g", xyzMM[i]);
    printf("\n");
  }
  {
    double tM[3][3], sM[3][3];
    for(int i=0; i<3; i++) for(int j=0; j<3; j++) tM[i][j]=rM[i][j];
    for(int i=0; i<3; i++) for(int j=0; j<3; j++) sM[i][j]=0.0;
    sM[0][0]=xyzMM[0]; sM[1][1]=xyzMM[1]; sM[2][2]=xyzMM[2];
    for(int i=0; i<3; i++)
      for(int j=0; j<3; j++)
        rM[i][j]= tM[i][0]*sM[0][j] + tM[i][1]*sM[1][j] + tM[i][2]*sM[2][j];
  }
 
  /* Fill most part of xform */
  for(int i=0; i<3; i++)
    for(int j=0; j<3; j++)
      xform[j+4*i]=rM[i][j];
  for(int i=0; i<3; i++)
    xform[3+4*i]=imgPos[i];
 
  /* Slice normalization vector not yet done */
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/
int dcmXformToQuatern(
  double *xform,
  double *quatern, 
  double *qoffset, 
  const int verbose
) {
  if(verbose>0) {printf("%s(xform[16], quatern[3], qoffset[3])\n", __func__); fflush(stdout);}
  if(xform==NULL) return(1);
  if(quatern!=NULL) for(int i=0; i<3; i++) quatern[i]=nan("");
  if(qoffset!=NULL) for(int i=0; i<3; i++) qoffset[i]=nan("");
 
  /* Copy array to matrix */
  double xformM[4][4];
  for(int i=0; i<4; i++)
    for(int j=0; j<4; j++)
      xformM[i][j]=xform[j+4*i];
  /* DICOM LPS to NIfTI RAS */
  for(int i=0; i<2; i++)
    for(int j=0; j<4; j++)
      xformM[i][j]=-xformM[i][j];
 
  /* Set quatern offsets */
  if(qoffset!=NULL) for(int i=0; i<3; i++) qoffset[i]=xformM[i][3];
 
  /* Ready, if quatern parameters not needed */
  if(quatern==NULL) return(0);
 
  // compute lengths
  double mt[3][3];
  for(int i=0; i<3; i++)
    for(int j=0; j<3; j++)
      mt[i][j]=xformM[i][j];
  double xd, yd, zd;
  xd=sqrt(mt[0][0]*mt[0][0] + mt[1][0]*mt[1][0] + mt[2][0]*mt[2][0]);
  yd=sqrt(mt[0][1]*mt[0][1] + mt[1][1]*mt[1][1] + mt[2][1]*mt[2][1]);
  zd=sqrt(mt[0][2]*mt[0][2] + mt[1][2]*mt[1][2] + mt[2][2]*mt[2][2]);
  if(!(xd>0.0)) {mt[0][0]=1.0; mt[1][0]=mt[2][0]=0.0; xd=1.0;}
  if(!(yd>0.0)) {mt[1][1]=1.0; mt[0][1]=mt[2][1]=0.0; yd=1.0;}
  if(!(zd>0.0)) {mt[2][2]=1.0; mt[0][2]=mt[1][2]=0.0; zd=1.0;}
  // normalize
  for(int i=0; i<3; i++) mt[i][0]/=xd;
  for(int i=0; i<3; i++) mt[i][1]/=yd;
  for(int i=0; i<3; i++) mt[i][2]/=zd;
  // check orthogonality (to be added)
 
  // determinant
  double dt;
  dt=  mt[0][0]*mt[1][1]*mt[2][2] - mt[0][0]*mt[2][1]*mt[1][2] - mt[1][0]*mt[0][1]*mt[2][2]
      +mt[1][0]*mt[2][1]*mt[0][2] + mt[2][0]*mt[0][1]*mt[1][2] - mt[2][0]*mt[1][1]*mt[0][2];
  if(verbose>3) printf("determinant=%g\n", dt);
  if(dt<0.0) for(int i=0; i<3; i++) mt[i][2]=-mt[i][2];
 
  // Calculate quatern parameters
  double a, b, c, d;
  a=1.0+mt[0][0]+mt[1][1]+mt[2][2];
  if(a>0.5) {
    a=0.5*sqrt(a);
    b=0.25*(mt[2][1]-mt[1][2])/a;
    c=0.25*(mt[0][2]-mt[2][0])/a;
    d=0.25*(mt[1][0]-mt[0][1])/a;
  } else {
    xd=1.0+mt[0][0]-(mt[1][1]+mt[2][2]);
    yd=1.0+mt[1][1]-(mt[0][0]+mt[2][2]);
    zd=1.0+mt[2][2]-(mt[0][0]+mt[1][1]);
    if(xd>1.0) {
      b=0.50*sqrt(xd);
      c=0.25*(mt[0][1]+mt[1][0])/b;
      d=0.25*(mt[0][2]+mt[2][0])/b;
      a=0.25*(mt[2][1]-mt[1][2])/b;
    } else if(yd>1.0) {
      c=0.50*sqrt(yd);
      b=0.25*(mt[0][1]+mt[1][0])/c;
      d=0.25*(mt[1][2]+mt[2][0])/c;
      a=0.25*(mt[0][2]-mt[2][0])/c;
    } else {
      d=0.50*sqrt(zd);
      b=0.25*(mt[0][2]+mt[2][0])/d;
      c=0.25*(mt[1][2]+mt[2][1])/d;
      a=0.25*(mt[1][0]-mt[0][1])/d;
    }
  }
  if(a<0.0) {b=-b; c=-c; d=-d;}
  quatern[0]=b;
  quatern[1]=c;
  quatern[2]=d;
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/