C++ ModelObject::CreateModelImage方法代码示例

//.........这里部分代码省略.........
    nColumns = options.nColumns;
    nRows = options.nRows;
  }
  nPixels_tot = nColumns * nRows;


  /* Read in PSF image, if supplied */
  if (options.psfImagePresent) {
    printf("Reading PSF image (\"%s\") ...\n", options.psfFileName.c_str());
    psfPixels = ReadImageAsVector(options.psfFileName, &nColumns_psf, &nRows_psf);
    if (psfPixels == NULL) {
      fprintf(stderr,  "\n*** ERROR: Unable to read PSF image file \"%s\"!\n\n", 
      			options.psfFileName.c_str());
      exit(-1);
    }
    nPixels_psf = nColumns_psf * nRows_psf;
    printf("naxis1 [# pixels/row] = %d, naxis2 [# pixels/col] = %d; nPixels_tot = %d\n", 
           nColumns_psf, nRows_psf, nPixels_psf);
  }
  else
    printf("* No PSF image supplied -- no image convolution will be done!\n");

  if (! options.subsamplingFlag)
    printf("* Pixel subsampling has been turned OFF.\n");

  /* Set up the model object */
  theModel = new ModelObject();
  theModel->SetDebugLevel(options.debugLevel);
  
  /* Add functions to the model object; also tells model object where function
     sets start */
  status = AddFunctions(theModel, functionList, functionBlockIndices, options.subsamplingFlag);
  if (status < 0) {
  	fprintf(stderr, "*** ERROR: Failure in AddFunctions!\n\n");
  	exit(-1);
  }

  
  // Add PSF image vector, if present (needs to be added prior to image data, so that
  // ModelObject can figure out proper internal model-image size
  if (options.psfImagePresent) {
    status = theModel->AddPSFVector(nPixels_psf, nColumns_psf, nRows_psf, psfPixels);
    if (status < 0) {
      fprintf(stderr, "*** ERROR: Failure in ModelObject::AddPSFVector!\n\n");
  	  exit(-1);
    }
  }

  /* Define the size of the requested model image */
  theModel->SetupModelImage(nColumns, nRows);
  theModel->PrintDescription();


  // Set up parameter vector(s), now that we know how many total parameters
  // there will be
  nParamsTot = theModel->GetNParams();
  printf("%d total parameters\n", nParamsTot);
  if (nParamsTot != (int)parameterList.size()) {
  	fprintf(stderr, "*** ERROR: number of input parameters (%d) does not equal", 
  	       (int)parameterList.size());
  	fprintf(stderr, " required number of parameters for specified functions (%d)!\n\n",
  	       nParamsTot);
  	exit(-1);
  }

  /* Copy parameters into C array and generate the model image */
  paramsVect = (double *) calloc(nParamsTot, sizeof(double));
  for (int i = 0; i < nParamsTot; i++)
    paramsVect[i] = parameterList[i];
  
  
  // Generate the image (including convolution, if requested), repeatedly
  gettimeofday(&timer_start, NULL);
  for (int ii = 0; ii < options.nIterations; ii++)
    theModel->CreateModelImage(paramsVect);
  gettimeofday(&timer_end, NULL);
  microsecs = timer_end.tv_usec - timer_start.tv_usec;
  time_elapsed = timer_end.tv_sec - timer_start.tv_sec + microsecs/1e6;
  printf("\nELAPSED TIME FOR %d ITERATIONS: %.6f sec\n", options.nIterations, time_elapsed);
  printf("Mean time per iteration = %.7f\n", time_elapsed/options.nIterations);

  

  /* Save model image: */
  printf("\nSaving output model image (\"%s\") ...\n", options.outputImageName.c_str());
  vector<string> commentStrings;
  SaveVectorAsImage(theModel->GetModelImageVector(), options.outputImageName, 
                    nColumns, nRows, commentStrings);

  printf("Done!\n\n");


  // Free up memory
  if (options.psfImagePresent)
    free(psfPixels);
  free(paramsVect);
  delete theModel;
  
  return 0;
}

//.........这里部分代码省略.........
  mp_par newParamLimit;
  for (int i = 0; i < nParamsTot; i++) {
    memset(&newParamLimit, 0, sizeof(mp_par));
    newParamLimit.fixed = paramLimits[i].fixed;
    if (paramLimits[i].fixed == 1) {
      printf("Fixed parameter detected (i = %d)\n", i);
      nFreeParams--;
    }
    newParamLimit.limited[0] = paramLimits[i].limited[0];
    newParamLimit.limited[1] = paramLimits[i].limited[1];
    newParamLimit.limits[0] = paramLimits[i].limits[0];
    newParamLimit.limits[1] = paramLimits[i].limits[1];
    parameterInfo.push_back(newParamLimit);
  }
  nDegFreedom = theModel->GetNValidPixels() - nFreeParams;
  printf("%d free parameters (%d degrees of freedom)\n", nFreeParams, nDegFreedom);

  // tell ModelObject about parameterInfo (mainly useful for printing-related methods)
  theModel->AddParameterInfo(parameterInfo);


  // OK, now we either print chi^2 value for the input parameters and quit, or
  // else call one of the solvers!
  if (options.printChiSquaredOnly) {
    printf("\n");
    fitStatus = 1;
    PrintFitStatistic(paramsVect, theModel, nFreeParams);
    printf("\n");
    // turn off saveing of parameter file
    options.saveBestFitParams = false;
  }
  else {
    // DO THE FIT!
    fitStatus = DispatchToSolver(options.solver, nParamsTot, nFreeParams, nStoredDataVals, 
    							paramsVect, parameterInfo, theModel, options.ftol, paramLimitsExist, 
    							options.verbose, &resultsFromSolver, options.nloptSolverName,
    							options.rngSeed);
    							
    PrintResults(paramsVect, theModel, nFreeParams, fitStatus, resultsFromSolver);
  }



  // OPTIONAL HANDLING OF BOOTSTRAP RESAMPLING HERE
  if ((options.doBootstrap) && (options.bootstrapIterations > 0)) {
    printf("\nNow doing bootstrap resampling (%d iterations) to estimate errors...\n",
           options.bootstrapIterations);
    printf("[NOT YET PROPERLY IMPLEMENTED!]\n");
    BootstrapErrors(paramsVect, parameterInfo, paramLimitsExist, theModel,
                    options.ftol, options.bootstrapIterations, nFreeParams);
  }
  
  
  
  if (options.saveBestFitParams) {
    printf("Saving best-fit parameters in file \"%s\"\n", options.outputParameterFileName.c_str());
    string  progNameVer = "profilefit ";
    progNameVer += VERSION_STRING;
    SaveParameters(paramsVect, theModel, options.outputParameterFileName,
                    programHeader, nFreeParams, options.solver, fitStatus, resultsFromSolver);
  }


  if (options.saveBestProfile) {
    theModel->CreateModelImage(paramsVect);
    double *modelProfile = (double *) calloc((size_t)nStoredDataVals, sizeof(double));
    int nPts = theModel->GetModelVector(modelProfile);
    if (nPts == nStoredDataVals) {
      printf("Saving model profile to %s...\n", options.modelOutputFileName.c_str());
      outputFile_ptr = fopen(options.modelOutputFileName.c_str(), "w");
      for (int i = 0; i < nPts; i++) {
        fprintf(outputFile_ptr, "\t%f\t%f\n", xVals[i], modelProfile[i]);
      }
      fclose(outputFile_ptr);
      printf("Done.\n");
    }
    else {
      printf("WARNING -- MISMATCH BETWEEN nStoredDataVals (main) and nDataVals (ModelObject1d)!\n");
      printf("NO PROFILE SAVED!\n");
    }
    free(modelProfile);
  }

  // Free up memory
  free(xVals);
  free(yVals);
  free(yWeights);
  if (maskAllocated)
    free(maskVals);
  if (options.psfPresent) {
    free(xVals_psf);
    free(yVals_psf);
  }
  free(paramsVect);
  free(paramErrs);
  delete theModel;
  
  printf("All done!\n\n");
  return 0;
}