imglhdv.c

Go to the documentation of this file.

/*****************************************************************************/
#include "tpcclibConfig.h"
/*****************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <string.h>
#include <math.h>
#include <time.h>
/*****************************************************************************/
#include "libtpccurveio.h"
#include "libtpcmodext.h"
#include "libtpcmisc.h"
#include "libtpcmodel.h"
#include "libtpcimgio.h"
#include "libtpcimgp.h"
/*****************************************************************************/
#define NNLS_N 4
#define MAX_MODEL_NR 2
enum {MODEL_SELECTION_WAIC, 
      MODEL_SELECTION_1, MODEL_SELECTION_2, MODEL_SELECTION_AIC};
/*****************************************************************************/
 
/*****************************************************************************/
static char *info[] = {
  "Computation of parametric image of distribution volume (DV) from dynamic",
  "PET image in ECAT, NIfTI, or Analyze format applying one- or two-tissue",
  "compartmental model with arterial plasma input.",
  "The compartmental models are transformed to general linear least squares",
  "functions, which are solved using Lawson-Hanson non-negative least squares",
  "(NNLS) algorithm (1). DV is estimated directly without division (2, 3, 4).",
  "Vascular volume is ignored here.",
  " ",
  "Dynamic PET image and plasma time-activity curve (PTAC) must be corrected",
  "for decay to the tracer injection time.",
  " ",
  "Usage: @P [Options] ptacfile imgfile dvfile",
  " ",
  "Options:",
  " -1 | -2 | -A | -0",
  "     With options -1 and -2 the one- or two-tissue compartment model can be",
  "     forced for all image voxels; otherwise both models are fitted, and",
  "     with option -A the model is selected based on AIC separately for each",
  "     voxel; by default (or option -0) Akaike weighted average of the model",
  "     parameters (5, 6) are reported.",
/*  " -K1=<filename>",
  "     Programs computes also a K1 image.",
*/
  " -m=<filename>",
  "     Programs writes the selected model number (1 or 2, or value between",
  "     1 and 2 as an image.",
  " -thr=<threshold%>",
  "     Pixels with AUC less than (threshold/100 x PTAC AUC) are set to zero",
  "     default is 0%",
  " -end=<Fit end time (min)>",
  "     Use data from 0 to end time; by default, model is fitted to all frames.",
  " -max=<Max value>",
  "     Upper limit for DV values.",
  " -stdoptions", // List standard options like --help, -v, etc
  " ",
  "The unit of voxel values in the DV image is (ml blood)/(ml tissue).",
  " ",
  "Example:",
  "  @P ua3818ap.kbq ua3818dy1.v ua3818dv.v",
  " ",
  "References:",
  "1. Lawson CL & Hanson RJ. Solving least squares problems.",
  "   Prentice-Hall, 1974, ISBN 0-89871-356-0.",
  "2. Zhou Y, Brasic J, Endres CJ, Kuwabara H, Kimes A, Contoreggi C, Maini A,",
  "   Ernst M, Wong DF. Binding potential image based statistical mapping for",
  "   detection of dopamine release by [11C]raclopride dynamic PET.",
  "   NeuroImage 2002;16(3):S91.",
  "3. Zhou Y, Brasic JR, Ye W, Dogan AS, Hilton J, Singer HS, Wong DF.",
  "   Quantification of cerebral serotonin binding in normal controls and",
  "   subjects with Tourette's syndrome using [11C]MDL 100,907 and",
  "   (+)[11C]McN 5652 dynamic PET with parametric imaging approach.",
  "   NeuroImage 2004;22(Suppl 2):T98.",
  "4. Hagelberg N, Aalto S, Kajander J, Oikonen V, Hinkka S, Någren K,",
  "   Hietala J, Scheinin H. Alfentanil increases cortical dopamine D2/D3",
  "   receptor binding in healthy subjects. Pain 2004;109:86-93.",
  "5. Turkheimer FE, Hinz R, Cunningham VJ. On the undecidability among",
  "   kinetic models: from model selection to model averaging. J Cereb Blood",
  "   Flow Metab 2003; 23: 490-498.",
  "6. Sederholm K. Model averaging with Akaike weights. TPCMOD0016 2003-04-07.",
  "   https://www.turkupetcentre.net/reports/tpcmod0016.pdf",
  " ",
  "See also: imgdv, imgbfbp, imgratio, img2tif, logan",
  " ",
  "Keywords: image, modelling, distribution volume, Vt, NNLS",
  0};
/*****************************************************************************/
 
/*****************************************************************************/
/* Turn on the globbing of the command line, since it is disabled by default in
   mingw-w64 (_dowildcard=0); in MinGW32 define _CRT_glob instead, if necessary;
   In Unix&Linux wildcard command line processing is enabled by default. */
/*
#undef _CRT_glob
#define _CRT_glob -1
*/
int _dowildcard = -1;
/*****************************************************************************/
 
/*****************************************************************************/
int main(int argc, char **argv)
{
  int      ai, help=0, version=0, verbose=1;
  char    *cptr;
  char     inpfile[FILENAME_MAX], petfile[FILENAME_MAX], dvfile[FILENAME_MAX];
  char     k1file[FILENAME_MAX], modfile[FILENAME_MAX];
  float    calcThreshold=0.0;
  double   upperLimit=-1.0;
  double   fittime=-1.0;
  int      model_selection=MODEL_SELECTION_WAIC;
  int      weight=0;
  int      ret;
  char     tmp[512];
 
 
  /*
   *  Get arguments
   */
  if(argc==1) {tpcPrintUsage(argv[0], info, stderr); return(1);}
  inpfile[0]=petfile[0]=dvfile[0]=k1file[0]=modfile[0]=(char)0;
  /* Get options */
  for(ai=1; ai<argc; ai++) if(*argv[ai]=='-') {
    if(tpcProcessStdOptions(argv[ai], &help, &version, &verbose)==0) continue;
    cptr=argv[ai]+1; if(*cptr=='-') cptr++; if(cptr==NULL) continue;
    if(strcasecmp(cptr, "A")==0) {
      model_selection=MODEL_SELECTION_AIC; continue;
    } else if(*cptr=='1') {
      model_selection=MODEL_SELECTION_1; continue;
    } else if(*cptr=='2') {
      model_selection=MODEL_SELECTION_2; continue;
    } else if(*cptr=='0') {
      model_selection=MODEL_SELECTION_WAIC; continue;
    } else if(strncasecmp(cptr, "K1=", 3)==0) {
      strlcpy(k1file, cptr+3, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "M=", 2)==0) {
      strlcpy(modfile, cptr+2, FILENAME_MAX); continue;
    } else if(strncasecmp(cptr, "THR=", 4)==0) {
      double v; ret=atof_with_check(cptr+4, &v);
      if(!ret && v>=0.0 && v<=200.0) {calcThreshold=(float)(0.01*v); continue;}
    } else if(strncasecmp(cptr, "END=", 4)==0) {
      fittime=atof_dpi(cptr+4); if(fittime>0.0) continue;
    } else if(*cptr=='w' || *cptr=='W') {
      weight=1; continue;
    } else if(strncasecmp(cptr, "MAX=", 4)==0) {
      upperLimit=atof_dpi(cptr+4); if(upperLimit>0.0) continue;
    }
    fprintf(stderr, "Error: invalid option '%s'.\n", argv[ai]);
    return(1);
  } else break;
 
  /* Print help or version? */
  if(help==2) {tpcHtmlUsage(argv[0], info, ""); return(0);}
  if(help) {tpcPrintUsage(argv[0], info, stdout); return(0);}
  if(version) {tpcPrintBuild(argv[0], stdout); return(0);}
  
  /* Process other arguments, starting from the first non-option */
  if(ai<argc) strlcpy(inpfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(petfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) strlcpy(dvfile, argv[ai++], FILENAME_MAX);
  if(ai<argc) {fprintf(stderr, "Error: invalid argument '%s'.\n", argv[ai]); return(1);}
  /* Did we get all the information that we need? */
  if(!dvfile[0]) {
    fprintf(stderr, "Error: missing command-line argument; use option --help\n");
    return(1);
  }
  if(strcasecmp(petfile, dvfile)==0 || strcasecmp(petfile, k1file)==0 ||
     strcasecmp(petfile, modfile)==0 || strcasecmp(inpfile, dvfile)==0)
  {
    fprintf(stderr, "Error: check the output filenames.\n");
    return(1);
  }
  /* In verbose mode print arguments and options */
  if(verbose>1) {
    printf("inpfile := %s\n", inpfile);
    printf("petfile := %s\n", petfile);
    printf("dvfile := %s\n", dvfile);
    if(k1file[0]) printf("k1file := %s\n", k1file);
    printf("calcThreshold :=%g\n", calcThreshold);
    if(upperLimit>0.0) printf("upperLimit := %g\n", upperLimit);
    printf("weight := %d\n", weight);
    if(fittime>0.0) printf("required_fittime := %g min\n", fittime);
    printf("model_selection := %d\n", model_selection);
    fflush(stdout);
  }
  if(verbose>8) IMG_TEST=verbose-8; else IMG_TEST=0;
 
 
  /*
   *  Read PET image and input TAC
   */
  if(verbose>0) printf("reading data files\n");
  DFT tac; dftInit(&tac); 
  IMG img; imgInit(&img);
  int dataNr;
  ret=imgReadModelingData(
    petfile, NULL, inpfile, NULL, NULL, &fittime, &dataNr, &img,
    NULL, &tac, 1, stdout, verbose-2, tmp);
  if(ret!=0) {
    fprintf(stderr, "Error: %s.\n", tmp);
    if(verbose>1) printf("  ret := %d\n", ret);
    return(2);
  }
  if(imgNaNs(&img, 1)>0)
    if(verbose>0) fprintf(stderr, "Warning: missing pixel values.\n");
  /* Set time unit to min, also for integrals in y2[] */
  if(tac.timeunit==TUNIT_SEC) {
    for(int fi=0; fi<tac.frameNr; fi++) tac.voi[0].y2[fi]/=60.0;
    for(int fi=0; fi<tac.frameNr; fi++) tac.voi[0].y3[fi]/=3600.0;
  }
  ret=dftTimeunitConversion(&tac, TUNIT_MIN);
  if(verbose>1) {
    printf("fittimeFinal := %g min\n", fittime);
    printf("dataNr := %d\n", dataNr);
  }
  /* Check that image is dynamic */
  if(dataNr<4) {
    fprintf(stderr, "Error: too few time frames for fitting.\n");
    if(verbose>0) imgInfo(&img);
    imgEmpty(&img); dftEmpty(&tac); return(2);
  }
 
 
  /* Add place for tissue TACs too */
  if(verbose>1) fprintf(stdout, "allocating working memory for pixel TACs\n");
  ret=dftAddmem(&tac, 1);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory.\n");
    if(verbose>0) printf("ret := %d\n", ret);
    imgEmpty(&img); dftEmpty(&tac); return(3);
  }
  strcpy(tac.voi[0].name, "input");
  strcpy(tac.voi[1].name, "tissue");
 
 
  /*
   *  Allocate result images, for DV, and for the selected model
   */
  if(verbose>1) fprintf(stdout, "allocating memory for parametric image data\n");
  IMG dvimg; imgInit(&dvimg);
  IMG mimg; imgInit(&mimg);
  ret=imgAllocateWithHeader(&dvimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(!ret) ret=imgAllocateWithHeader(&mimg, img.dimz, img.dimy, img.dimx, 1, &img);
  if(ret) {
    fprintf(stderr, "Error: cannot allocate memory for result image.\n");
    imgEmpty(&img); dftEmpty(&tac); imgEmpty(&dvimg); imgEmpty(&mimg); return(4);
  }
  /* and set the more or less necessary things */
  dvimg.start[0]=mimg.start[0]=0.0; 
  dvimg.end[0]=mimg.end[0]=60.*fittime;
  dvimg.unit=CUNIT_ML_PER_ML; /* mL/mL */
  mimg.unit=CUNIT_UNITLESS;
 
 
  /*
   *  Allocate memory required by NNLS
   */
  if(verbose>1) fprintf(stdout, "allocating memory for NNLS\n");
  int      n, m, nnls_n, nnls_m, nnls_index[NNLS_N];
  double  *nnls_a[NNLS_N], *nnls_b, *nnls_zz, nnls_x[NNLS_N], *nnls_mat,
           nnls_wp[NNLS_N], *dptr, nnls_rnorm;
  nnls_n=NNLS_N; nnls_m=dataNr;
  nnls_mat=(double*)malloc(((nnls_n+2)*nnls_m)*sizeof(double));
  if(nnls_mat==NULL) {
    fprintf(stderr, "Error: cannot allocate memory for NNLS.\n");
    imgEmpty(&img); dftEmpty(&tac); imgEmpty(&dvimg); imgEmpty(&mimg);
    return(4);
  }
  for(n=0, dptr=nnls_mat; n<nnls_n; n++) {nnls_a[n]=dptr; dptr+=nnls_m;}
  nnls_b=dptr; dptr+=nnls_m; nnls_zz=dptr;
 
  /* Copy weights if available */
  /* or set them to frame lengths */
  if(weight==1 && img.isWeight==0) {
    for(m=0; m<nnls_m; m++) img.weight[m]=img.end[m]-img.start[m];
    img.isWeight=1;
  }
  /* Compute weights for NNLS */
  if(img.isWeight) {
    for(m=0; m<nnls_m; m++) tac.w[m]=img.weight[m];
    for(m=0; m<nnls_m; m++) {
      if(tac.w[m]<=1.0e-20) tac.w[m]=0.0; else tac.w[m]=sqrt(tac.w[m]);
    }
  }
 
 
  /*
   *  Compute pixel-by-pixel
   */
  /* Determine the threshold */
  double threshold=calcThreshold*tac.voi[0].y2[dataNr-1];
  if(verbose>1) printf("threshold_AUC := %g\n", threshold);
  if(verbose>0) fprintf(stdout, "computing NNLS pixel-by-pixel\n");
  double f, g;
  double aic[MAX_MODEL_NR], w[MAX_MODEL_NR], ss[MAX_MODEL_NR], dv[MAX_MODEL_NR];
  int cm2ok=0;
  for(int pi=0; pi<img.dimz; pi++) {
    if(verbose>1) printf("computing plane %d\n", img.planeNumber[pi]);
    else if(verbose>0) fprintf(stdout, ".");
    if(verbose>0) fflush(stdout);
    for(int yi=0; yi<img.dimy; yi++) {
      for(int xi=0; xi<img.dimx; xi++) {
        /* Set pixel results to zero */
        dvimg.m[pi][yi][xi][0]=0.0; 
        mimg.m[pi][yi][xi][0]=0.0;
        /* Copy and integrate pixel curve */
        for(int fi=0; fi<tac.frameNr; fi++)
          tac.voi[1].y[fi]=img.m[pi][yi][xi][fi];
        ret=petintegral(tac.x1, tac.x2, tac.voi[1].y, tac.frameNr, tac.voi[1].y2, tac.voi[1].y3);
        if(ret) continue;
        /* if AUC at the end is less than threshold value, then do nothing */
        if(tac.voi[1].y2[dataNr-1]<threshold) continue;
 
        /* Fit two-tissue model, if required */
        if(model_selection==MODEL_SELECTION_AIC || 
           model_selection==MODEL_SELECTION_2 || 
           model_selection==MODEL_SELECTION_WAIC)
        {
          nnls_n=4;
          /* Fill the NNLS data matrix */
          for(m=0; m<nnls_m; m++) {
            nnls_a[0][m]=tac.voi[0].y3[m];
            nnls_a[1][m]=-tac.voi[1].y2[m];
            nnls_a[2][m]=tac.voi[0].y2[m];
            nnls_a[3][m]=-tac.voi[1].y[m];
            nnls_b[m]=tac.voi[1].y3[m];
          }
          if(img.isWeight) { 
            for(m=0; m<nnls_m; m++) {
              nnls_b[m]*=tac.w[m];
              for(n=0; n<nnls_n; n++) nnls_a[n][m]*=tac.w[m];
            }
          }
          /* NNLS */
          ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm,
                   nnls_wp, nnls_zz, nnls_index);
          if(ret>1) { /* no solution is possible */
            for(n=0; n<nnls_n; n++) nnls_x[n]=0.0;
            nnls_rnorm=0.0;
          } else if(ret==1) { /* max iteration count exceeded */ }
          /* Get DV as the 1st linear coefficient */
          dv[1]=nnls_x[0];
          /* Calculate Sum-of-Squares */
          for(m=0; m<nnls_m; m++) {
            nnls_a[0][m]=tac.voi[0].y3[m];
            nnls_a[1][m]=-tac.voi[1].y2[m];
            nnls_a[2][m]=tac.voi[0].y2[m];
            nnls_a[3][m]=-tac.voi[1].y[m];
            nnls_b[m]=tac.voi[1].y3[m];
          }
          if(img.isWeight) {
            for(m=0; m<nnls_m; m++) {
              nnls_b[m]*=tac.w[m];
              for(n=0; n<nnls_n; n++) nnls_a[n][m]*=tac.w[m];
            }
          }
          for(m=0, ss[1]=0.0; m<nnls_m; m++) {
            for(n=0, f=0.0; n<nnls_n; n++) f+=nnls_x[n]*nnls_a[n][m];
            g=f; f-=nnls_b[m]; ss[1]+=f*f;
            if(verbose>17 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3)
              printf("%10.3f :  %e  %e - %e = %g\n",
                     tac.x[m], tac.voi[1].y[m], tac.voi[1].y3[m], g, f);
          }
          /* Calculate AIC */
          aic[1]=aicSS(ss[1], nnls_m, nnls_n);
          if(verbose>15 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
            for(n=0; n<nnls_n; n++) printf(" p[%d]=%g", n, nnls_x[n]);
            printf(" -> SS=%g AIC=%g\n", ss[1], aic[1]);
          }
          /* Compute how many of the parameters are not zeros */
          for(n=0, m=nnls_n; n<nnls_n; n++) if(nnls_x[n]<1.0e-10) m--;
          /* Check if results were ok */
          if(m==nnls_n) cm2ok=1; else cm2ok=0;
        } /* done with two-tissue model */
 
        /* Fit one-tissue model, if required */
        if(model_selection==MODEL_SELECTION_AIC ||
           model_selection==MODEL_SELECTION_1 ||
           model_selection==MODEL_SELECTION_WAIC)
        {
          nnls_n=2;
          /* Fill the NNLS data matrix */
          for(m=0; m<nnls_m; m++) {
            nnls_a[0][m]=tac.voi[0].y3[m];
            nnls_a[1][m]=-tac.voi[1].y2[m];
            nnls_b[m]=tac.voi[1].y3[m];
          }
          if(img.isWeight) {
            for(m=0; m<nnls_m; m++) {
              nnls_b[m]*=tac.w[m];
              for(n=0; n<nnls_n; n++) nnls_a[n][m]*=tac.w[m];
            }
          }
          /* NNLS */
          ret=nnls(nnls_a, nnls_m, nnls_n, nnls_b, nnls_x, &nnls_rnorm,
                   nnls_wp, nnls_zz, nnls_index);
          if(ret>1) { /* no solution is possible */
            for(n=0; n<nnls_n; n++) nnls_x[n]=0.0;
            nnls_rnorm=0.0;
          } else if(ret==1) { /* max iteration count exceeded */ }
          /* Get DV as the 1st linear coefficient */
          dv[0]=nnls_x[0];
          /* Calculate Sum-of-Squares */
          for(m=0; m<nnls_m; m++) {
            nnls_a[0][m]=tac.voi[0].y3[m];
            nnls_a[1][m]=-tac.voi[1].y2[m];
            nnls_b[m]=tac.voi[1].y3[m];
          }
          if(img.isWeight) {
            for(m=0; m<nnls_m; m++) {
              nnls_b[m]*=tac.w[m];
              for(n=0; n<nnls_n; n++) nnls_a[n][m]*=tac.w[m];
            }
          }
          for(m=0, ss[0]=0.0; m<nnls_m; m++) {
            for(n=0, f=0.0; n<nnls_n; n++) f+=nnls_x[n]*nnls_a[n][m];
            g=f; f-=nnls_b[m]; ss[0]+=f*f;
            if(verbose>17 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3)
              printf("%10.3f :  %e  %e - %e = %g\n",
                     tac.x[m], tac.voi[1].y[m], tac.voi[1].y3[m], g, f);
          }
          /* Calculate AIC */
          aic[0]=aicSS(ss[0], nnls_m, nnls_n);
          if(verbose>15 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3) {
            for(n=0; n<nnls_n; n++) printf(" p[%d]=%g", n, nnls_x[n]);
            printf(" -> SS=%g AIC=%g\n", ss[0], aic[0]);
          }
        } /* done with one-tissue model */
 
        /* Select the model, or calculate average */
        switch(model_selection) {
        case MODEL_SELECTION_WAIC:
          ret=aicWeights(aic, w, 2);
          if(verbose>14 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3)
            printf(" Model weights are %g %g\n", w[0], w[1]);
          /* Compute weighted DV and other 'parameters' */
          dvimg.m[pi][yi][xi][0]=aicWeightedAvg(w, dv, 2);
          mimg.m[pi][yi][xi][0]=(float)(1.0*w[0]+2.0*w[1]);
          if(verbose>13 && pi==img.dimz/2 && yi==img.dimy/3 && xi==img.dimx/3)
            printf(" DV1=%g DV2=%g -> Weighted DV=%g\n", dv[0], dv[1], dvimg.m[pi][yi][xi][0]);
          break;
        case MODEL_SELECTION_1:
          dvimg.m[pi][yi][xi][0]=dv[0];
          mimg.m[pi][yi][xi][0]=(float)1.0;
          break;
        case MODEL_SELECTION_2:
          dvimg.m[pi][yi][xi][0]=dv[1];
          mimg.m[pi][yi][xi][0]=(float)2.0;
          break;
        case MODEL_SELECTION_AIC:
          if(cm2ok && aic[1]<aic[0]) {
            dvimg.m[pi][yi][xi][0]=dv[1];
            mimg.m[pi][yi][xi][0]=2.0;
          } else {
            dvimg.m[pi][yi][xi][0]=dv[0];
            mimg.m[pi][yi][xi][0]=1.0;
          }
          break;
        }
 
      } /* next column */
    } /* next row */
  } /* next plane */
  if(verbose>0) fprintf(stdout, " done.\n");
  /* No need for NNLS data anymore */
  free(nnls_mat);
  /* No need for dynamic image or input data anymore */
  imgEmpty(&img);
  dftEmpty(&tac);
 
 
  /*
   *  Filter the parametric images if required
   */
  /* Remove values that are higher than the user-defined limit */
  if(upperLimit>0.0) {
    if(verbose>0) fprintf(stdout, "filtering out pixel values exceeding %g\n", upperLimit);
    imgCutoff(&dvimg, (float)upperLimit, 0);
  }
 
 
  /*
   *  Save parametric image
   */
  ret=backupExistingFile(dvfile, NULL, tmp); if(ret!=0) {
    fprintf(stderr, "Error: %s\n", tmp); 
    imgEmpty(&dvimg); imgEmpty(&mimg);
    return(11);
  }
  if(verbose>1) printf("%s\n", tmp);
  ret=imgWrite(dvfile, &dvimg);
  if(ret) {
    fprintf(stderr, "Error: %s\n", dvimg.statmsg);
    imgEmpty(&dvimg); imgEmpty(&mimg);
    return(11);
  }
  if(verbose>0) fprintf(stdout, "DV image %s saved.\n", dvfile);
 
  /*
   *  Save selected model as image, if requested
   */
  if(modfile[0]) {
    ret=backupExistingFile(modfile, NULL, tmp); if(ret!=0) {
      fprintf(stderr, "Error: %s\n", tmp); 
      imgEmpty(&dvimg); imgEmpty(&mimg);
      return(13);
    }
    if(verbose>1) printf("%s\n", tmp);
    ret=imgWrite(modfile, &mimg);
    if(ret) {
      fprintf(stderr, "Error: %s\n", mimg.statmsg);
      imgEmpty(&dvimg); imgEmpty(&mimg);
      return(13);
    }
    if(verbose>0) fprintf(stdout, "Model selection image %s saved.\n", modfile);
  }
 
 
  imgEmpty(&dvimg); imgEmpty(&mimg);
 
  return(0);
}
/*****************************************************************************/
 
/*****************************************************************************/