C++ Projection::ToWorldX方法代码示例

  /**
   * Set the output cube to specified file name and specified input images
   * and output attributes and lat,lons
   */
  Isis::Cube *ProcessMapMosaic::SetOutputCube(const QString &inputFile,
      double xmin, double xmax, double ymin, double ymax,
      double slat, double elat, double slon, double elon, int nbands,
      CubeAttributeOutput &oAtt, const QString &mosaicFile) {
    Pvl fileLab(inputFile);
    PvlGroup &mapping = fileLab.findGroup("Mapping", Pvl::Traverse);

    mapping["UpperLeftCornerX"] = toString(xmin);
    mapping["UpperLeftCornerY"] = toString(ymax);
    mapping.addKeyword(PvlKeyword("MinimumLatitude", toString(slat)), Pvl::Replace);
    mapping.addKeyword(PvlKeyword("MaximumLatitude", toString(elat)), Pvl::Replace);
    mapping.addKeyword(PvlKeyword("MinimumLongitude", toString(slon)), Pvl::Replace);
    mapping.addKeyword(PvlKeyword("MaximumLongitude", toString(elon)), Pvl::Replace);

    Projection *firstProj = ProjectionFactory::CreateFromCube(fileLab);
    int samps = (int)(ceil(firstProj->ToWorldX(xmax) - firstProj->ToWorldX(xmin)) + 0.5);
    int lines = (int)(ceil(firstProj->ToWorldY(ymin) - firstProj->ToWorldY(ymax)) + 0.5);
    delete firstProj;

    if (p_createMosaic) {
      Pvl newMap;
      newMap.addGroup(mapping);

      // Initialize the mosaic
      CubeAttributeInput inAtt;

      ProcessByLine p;
      p.SetInputCube(inputFile, inAtt);
      p.PropagateHistory(false);
      p.PropagateLabels(false);
      p.PropagateTables(false);
      p.PropagatePolygons(false);
      p.PropagateOriginalLabel(false);

      // If track set, create the origin band
      if (GetTrackFlag()) {
        nbands += 1;
      }
      // For average priority, get the new band count
      else if (GetImageOverlay() == AverageImageWithMosaic) {
        nbands *= 2;
      }

      Cube *ocube = p.SetOutputCube(mosaicFile, oAtt, samps, lines, nbands);
      p.Progress()->SetText("Initializing mosaic");
      p.ClearInputCubes();
      p.StartProcess(ProcessMapMosaic::FillNull);

      // CreateForCube created some keywords in the mapping group that needs to be added
      ocube->putGroup(newMap.findGroup("Mapping", Pvl::Traverse));
      p.EndProcess();
    }

    Cube *mosaicCube = new Cube();
    mosaicCube->open(mosaicFile, "rw");
    mosaicCube->addCachingAlgorithm(new UniqueIOCachingAlgorithm(2));

    AddOutputCube(mosaicCube);
    return mosaicCube;
  }

  /**
   * Mosaic Processing method, returns false if the cube is not inside the mosaic
   */
  bool ProcessMapMosaic::StartProcess(QString inputFile) {
    if (InputCubes.size() != 0) {
      QString msg = "Input cubes already exist; do not call SetInputCube when using ";
      msg += "ProcessMosaic::StartProcess(QString)";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    if (OutputCubes.size() == 0) {
      QString msg = "An output cube must be set before calling StartProcess";
      throw IException(IException::Programmer, msg, _FILEINFO_);
    }

    CubeAttributeInput inAtt(inputFile);
    Cube *inCube = ProcessMosaic::SetInputCube(inputFile, inAtt);

    Cube *mosaicCube = OutputCubes[0];
    Projection *iproj = inCube->projection();
    Projection *oproj = mosaicCube->projection();
    int nsMosaic = mosaicCube->sampleCount();
    int nlMosaic = mosaicCube->lineCount();

    if (*iproj != *oproj) {
      QString msg = "Mapping groups do not match between cube [" + inputFile + "] and mosaic";
      throw IException(IException::User, msg, _FILEINFO_);
    }

    int outSample, outSampleEnd, outLine, outLineEnd;
    outSample = (int)(oproj->ToWorldX(iproj->ToProjectionX(1.0)) + 0.5);
    outLine   = (int)(oproj->ToWorldY(iproj->ToProjectionY(1.0)) + 0.5);

    int ins = InputCubes[0]->sampleCount();
    int inl =  InputCubes[0]->lineCount();
    outSampleEnd = outSample + ins;
    outLineEnd   = outLine + inl;

    bool wrapPossible = iproj->IsEquatorialCylindrical();
    int worldSize = 0;
    if (wrapPossible) {
      // Figure out how many samples 360 degrees is
      wrapPossible = wrapPossible && oproj->SetUniversalGround(0, 0);
      int worldStart = (int)(oproj->WorldX() + 0.5);
      wrapPossible = wrapPossible && oproj->SetUniversalGround(0, 180);
      int worldEnd = (int)(oproj->WorldX() + 0.5);

      worldSize = abs(worldEnd - worldStart) * 2;

      wrapPossible = wrapPossible && (worldSize > 0);

      // This is EquatorialCylindrical, so shift to the left all the way
      if (wrapPossible) {
        // While some data would still be put in the mosaic, move left
        //  >1 for end because 0 still means no data, whereas 1 means 1 line of data
        while (outSampleEnd - worldSize > 1) {
          outSample -= worldSize;
          outSampleEnd -= worldSize;
        }
        // Now we have the sample range to the furthest left
      }
    }

    // Check overlaps of input image along the mosaic edges before
    // calling ProcessMosaic::StartProcess
    // Left edge
    if (outSample < 1) {
      ins = ins + outSample - 1;
    }

    // Top edge
    if (outLine < 1) {
      inl = inl + outLine - 1;
    }

    // Right edge
    if ((outSample + ins - 1) > nsMosaic) {
      ins = nsMosaic - outSample + 1;
    }

    // Bottom edge
    if ((outLine + inl - 1) > nlMosaic) {
      inl = nlMosaic - outLine + 1;
    }

    if (outSampleEnd < 1 || outLineEnd < 1 || outSample > nsMosaic || outLine > nlMosaic || ins < 1 || inl < 1) {
      // Add a PvlKeyword naming which files are not included in output mosaic
      ClearInputCubes();
      return false;
    }
    else {
      // Place the input in the mosaic
      Progress()->SetText("Mosaicking " + FileName(inputFile).name());

      try {
        do {
          int outBand = 1;
          ProcessMosaic::StartProcess(outSample, outLine, outBand);

          // Increment for projections where occurrances may happen multiple times
//.........这里部分代码省略.........

//.........这里部分代码省略.........
    if (runXY) {
      double startX = DBL_MAX;
      double endX = DBL_MIN;
      double startY = DBL_MAX;
      double endY =  DBL_MIN;
      for (int i=0; i<(int)clist.size(); i++) {
        Projection *proj = clist[i]->Projection();
        proj->SetWorld(0.5,0.5);
        if (i==0) {
          startX = proj->XCoord();
          endY = proj->YCoord();
        }
        else {
          if (proj->XCoord() < startX) startX =  proj->XCoord();
          if (proj->YCoord() > endY) endY = proj->YCoord();
        }
        Pvl *p = clist[i]->Label();
        double nlines = p->FindGroup("Dimensions",Pvl::Traverse)["Lines"];
        double nsamps = p->FindGroup("Dimensions",Pvl::Traverse)["Samples"];

        proj->SetWorld((nsamps+0.5),(nlines+0.5));
        if (i==0) {
          endX = proj->XCoord();
          startY = proj->YCoord();
        }
        else {
          if (proj->XCoord() > endX) endX =  proj->XCoord();
          if (proj->YCoord() < startY) startY = proj->YCoord();
        }
      }

      double avgX = (startX + endX) / 2;
      double avgY = (startY + endY) / 2;
      double sample = proj->ToWorldX(avgX);
      double line = proj->ToWorldY(avgY);

      for (int i=0; i<(int)clist.size(); i++) {
        Camera *cam = clist[i]->Camera();
        if (cam->SetImage(sample,line)) {
          Cemiss = cam->EmissionAngle();
          Cphase = cam->PhaseAngle();
          Cincid = cam->IncidenceAngle();
          ClocalSolTime = cam->LocalSolarTime();
          CsolarLong = cam->SolarLongitude();
          CsunAzimuth = cam->SunAzimuth();
          CnorthAzimuth = cam->NorthAzimuth();
          runXY = false;
          break;
        }
      }
    }
    if (runXY) {
      string msg = "Camera did not intersect images to gather stats";
      throw iException::Message(iException::User,msg,_FILEINFO_);
    }

    // get the min and max SCLK values ( do this with string comp.)
    // get the value from the original label blob
    string startClock;
    string stopClock;
    string startTime;
    string stopTime;
    for (int i=0; i<(int)clist.size(); i++) {
      OriginalLabel origLab;
      clist[i]->Read(origLab);
      PvlGroup timegrp = origLab.ReturnLabels().FindGroup("TIME_PARAMETERS",Pvl::Traverse);

/**
 * Constructs an OverlapStatistics object.  Compares the two input cubes and
 * finds where they overlap.
 *
 * @param x The first input cube
 * @param y The second input cube
 * @param progressMsg (Default value of "Gathering Overlap Statistics") Text
 *         for indicating progress during statistic gathering
 * @param sampPercent (Default value of 100.0) Sampling percent, or the percentage
 *       of lines to consider during the statistic gathering procedure
 *
 * @throws Isis::iException::User - All images must have the same number of
 *                                  bands
 */
OverlapStatistics::OverlapStatistics(Isis::Cube &x, Isis::Cube &y,
                                     std::string progressMsg, double sampPercent) {
    // Test to ensure sampling percent in bound
    if (sampPercent <= 0.0 || sampPercent > 100.0) {
        string msg = "The sampling percent must be a decimal (0.0, 100.0]";
        throw iException::Message(iException::Programmer,msg,_FILEINFO_);
    }

    p_sampPercent = sampPercent;

    // Extract filenames and band number from cubes
    p_xFile = x.Filename();
    p_yFile = y.Filename();

    // Make sure number of bands match
    if (x.Bands() != y.Bands()) {
        string msg = "Number of bands do not match between cubes [" +
                     p_xFile.Name() + "] and [" + p_yFile.Name() + "]";
        throw iException::Message(iException::User,msg,_FILEINFO_);
    }
    p_bands = x.Bands();
    p_stats.resize(p_bands);

    //Create projection from each cube
    Projection *projX = x.Projection();
    Projection *projY = y.Projection();

    // Test to make sure projection parameters match
    if (*projX != *projY) {
        string msg = "Mapping groups do not match between cubes [" +
                     p_xFile.Name() + "] and [" + p_yFile.Name() + "]";
        throw iException::Message(iException::Programmer,msg,_FILEINFO_);
    }

    // Figure out the x/y range for both images to find the overlap
    double Xmin1 = projX->ToProjectionX(0.5);
    double Ymax1 = projX->ToProjectionY(0.5);
    double Xmax1 = projX->ToProjectionX(x.Samples()+0.5);
    double Ymin1 = projX->ToProjectionY(x.Lines()+0.5);

    double Xmin2 = projY->ToProjectionX(0.5);
    double Ymax2 = projY->ToProjectionY(0.5);
    double Xmax2 = projY->ToProjectionX(y.Samples()+0.5);
    double Ymin2 = projY->ToProjectionY(y.Lines()+0.5);

    // Find overlap
    if ((Xmin1<Xmax2) && (Xmax1>Xmin2) && (Ymin1<Ymax2) && (Ymax1>Ymin2)) {
        double minX = Xmin1 > Xmin2 ? Xmin1 : Xmin2;
        double minY = Ymin1 > Ymin2 ? Ymin1 : Ymin2;
        double maxX = Xmax1 < Xmax2 ? Xmax1 : Xmax2;
        double maxY = Ymax1 < Ymax2 ? Ymax1 : Ymax2;

        // Find Sample range of the overlap
        p_minSampX = (int)(projX->ToWorldX(minX) + 0.5);
        p_maxSampX = (int)(projX->ToWorldX(maxX) + 0.5);
        p_minSampY = (int)(projY->ToWorldX(minX) + 0.5);
        p_maxSampY = (int)(projY->ToWorldX(maxX) + 0.5);
        p_sampRange = p_maxSampX - p_minSampX + 1;

        // Test to see if there was only sub-pixel overlap
        if (p_sampRange <= 0) return;

        // Find Line range of overlap
        p_minLineX = (int)(projX->ToWorldY(maxY) + 0.5);
        p_maxLineX = (int)(projX->ToWorldY(minY) + 0.5);
        p_minLineY = (int)(projY->ToWorldY(maxY) + 0.5);
        p_maxLineY = (int)(projY->ToWorldY(minY) + 0.5);
        p_lineRange = p_maxLineX - p_minLineX + 1;

        // Print percent processed
        Progress progress;
        progress.SetText(progressMsg);

        int linc = (int)(100.0 / sampPercent + 0.5); // Calculate our line increment

        // Define the maximum number of steps to be our line range divided by the
        // line increment, but if they do not divide evenly, then because of
        // rounding, we need to do an additional step for each band
        int maxSteps = (int)(p_lineRange / linc + 0.5);

        if (p_lineRange % linc != 0) maxSteps += 1;
        maxSteps *= p_bands;


        progress.SetMaximumSteps(maxSteps);
        progress.CheckStatus();
//.........这里部分代码省略.........