C++ Image::At方法代码示例

TEST(MergeAcq_Test, TopMergeTest) {
  Image top;
  Image bottom;
  Image combo;
  MergeAcq merger;
  // Load up our test files
  const char *file = "test.dat";
  ION_ASSERT(top.LoadRaw(file), "Couldn't load file.");
  ION_ASSERT(bottom.LoadRaw(file), "Couldn't load file.");
  merger.SetFirstImage(&bottom);
  merger.SetSecondImage(&top, bottom.GetRows(), 0); // starting vertically raised but columns the same.
  // Merge the two images into a single image
  merger.Merge(combo);
  cout << "Done." << endl;
  // Write out the merged image for some testing
  Acq acq;
  acq.SetData(&combo);
  acq.WriteVFC("combo.dat", 0, 0, combo.GetCols(), combo.GetRows());
  Image test;
  ION_ASSERT(test.LoadRaw("combo.dat"), "Couldn't load file.");
  // Check to make sure that the combined file has same values as original
  for (int row = 0; row < top.GetRows(); row++) {
    for (int col = 0; col < top.GetCols(); col++) {
      for (int frame = 0; frame < top.GetFrames(); frame++) {
        EXPECT_EQ(top.At(row, col, frame), bottom.At(row,col,frame));
        short orig = top.At(row, col, frame);
        short combined = test.At(row+top.GetRows(),col,frame);
        EXPECT_EQ(orig, combined);
        combined = test.At(row,col,frame);
        EXPECT_EQ(orig, combined);
      }
    }
  }
}

void BFReference::CalcSignalReference2(const std::string &datFile, const std::string &bgFile,
				      Mask &mask, int traceFrame) {
  Image bfImg;
  Image bfBkgImg;
  bfImg.SetImgLoadImmediate (false);
  bfBkgImg.SetImgLoadImmediate (false);
  bool loaded = bfImg.LoadRaw(datFile.c_str());
  bool bgLoaded = bfBkgImg.LoadRaw(bgFile.c_str());
  if (!loaded) {
    ION_ABORT("*Error* - No beadfind file found, did beadfind run? are files transferred?  (" + datFile + ")");
  }
  if (!bgLoaded) {
    ION_ABORT("*Error* - No beadfind background file found, did beadfind run? are files transferred?  (" + bgFile + ")");
  }
  const RawImage *raw = bfImg.GetImage();
  
  assert(raw->cols == GetNumCol());
  assert(raw->rows == GetNumRow());
  assert(raw->cols == mask.W());
  assert(raw->rows == mask.H());
  int StartFrame = bfImg.GetFrame(-663); //5
  int EndFrame = bfImg.GetFrame(350); //20
  int NNinnerx = 1, NNinnery = 1, NNouterx = 12, NNoutery = 8;
  cout << "DC start frame: " << StartFrame << " end frame: " << EndFrame << endl;
  bfImg.FilterForPinned(&mask, MaskEmpty, false);
  bfImg.XTChannelCorrect();
  // bfImg.XTChannelCorrect(&mask);
  Traces trace;  
  trace.Init(&bfImg, &mask, FRAMEZERO, FRAMELAST, FIRSTDCFRAME,LASTDCFRAME);
  bfImg.Normalize(StartFrame, EndFrame);
  if (mDoRegionalBgSub) {
     trace.SetMeshDist(0);
  }
  trace.CalcT0(true);
  if (mDoRegionalBgSub) {
    GridMesh<float> grid;
    grid.Init(raw->rows, raw->cols, mRegionYSize, mRegionXSize);
    int numBin = grid.GetNumBin();
    int rowStart = -1, rowEnd = -1, colStart = -1, colEnd = -1;
    for (int binIx = 0; binIx < numBin; binIx++) {
      cout << "BG Subtract Region: " << binIx << endl;
      grid.GetBinCoords(binIx, rowStart, rowEnd, colStart, colEnd);
      Region reg;
      reg.row = rowStart;
      reg.h = rowEnd - rowStart;
      reg.col = colStart;
      reg.w = colEnd - colStart;
      bfImg.BackgroundCorrectRegion(&mask, reg, MaskAll, MaskEmpty, NNinnerx, NNinnery, NNouterx, NNoutery, NULL);
    }
  }
  else {
    bfImg.BackgroundCorrect(&mask, MaskEmpty, MaskEmpty, NNinnerx, NNinnery, NNouterx, NNoutery, NULL);
  }
  int length = GetNumRow() * GetNumCol();
  mBfMetric.resize(length, std::numeric_limits<double>::signaling_NaN());
  for (int wIx = 0; wIx < length; wIx++) {
    if (mask[wIx] & MaskExclude || mask[wIx] & MaskPinned) 
      continue;
    int t0 = (int)trace.GetT0(wIx);
    mBfMetric[wIx] = 0;
    float zSum  = 0;
    int count = 0;
    for (int fIx = min(t0-20, 0); fIx < t0-10; fIx++) {
      zSum += bfImg.At(wIx,fIx);
      count ++;
    }
    for (int fIx = t0+3; fIx < t0+15; fIx++) {
      mBfMetric[wIx] += (bfImg.At(wIx,fIx) - (zSum / count));
    }
  }
  bfImg.Close();
  for (int i = 0; i < length; i++) {
    if (mask[i] & MaskExclude || mWells[i] == Exclude) {
      mWells[i] = Exclude;
    }
    else {
      mask[i] = MaskIgnore;
    }
  }
  cout << "Filling reference. " << endl;
  FillInReference(mWells, mBfMetric, mGrid, mMinQuantile, mMaxQuantile, mNumEmptiesPerRegion);
  for (int i = 0; i < length; i++) {
    if (mWells[i] == Reference) {
      mask[i] = MaskEmpty;
    }
  }
}

void BFReference::FilterRegionOutliers(Image &bfImg, Mask &mask, float iqrThreshold, 
                                       int rowStart, int rowEnd, int colStart, int colEnd) {

  const RawImage *raw = bfImg.GetImage();
  /* Figure out how many wells are not pinned/excluded right out of the gate. */
  int okCount = 0;
  for (int r = rowStart; r < rowEnd; r++) {
    for (int c = colStart; c < colEnd; c++) {
      int idx = r * raw->cols + c;
      if (!mask.Match(c, r, MaskPinned) && !mask.Match(c,r,MaskExclude) && 
          mWells[idx] != Exclude && mWells[idx] != Filtered) {
        okCount++;
      }
    }
  }
  /* If not enough, just mark them all as bad. */
  if (okCount <= MIN_OK_WELLS || (mNumEmptiesPerRegion > 0 && okCount < mNumEmptiesPerRegion)) {
    for (int r = rowStart; r < rowEnd; r++) {
      for (int c = colStart; c < colEnd; c++) {
        mWells[r * raw->cols + c] = Exclude;
      }
    }
    return;
  }

  // Make a mtrix for our region
  mTraces.set_size(okCount, raw->frames ); // wells in row major order by frames
  int count = 0;

  vector<int> mapping(mTraces.n_rows, -1);
  for (int r = rowStart; r < rowEnd; r++) {
    for (int c = colStart; c < colEnd; c++) {
      int idx = r * raw->cols + c;
      if (!mask.Match(c, r, MaskPinned) && !mask.Match(c,r,MaskExclude) && 
          mWells[idx] != Exclude && mWells[idx] != Filtered) {
        for (int f = 0; f < raw->frames; f++) {
          mTraces.at(count,f) = bfImg.At(r,c,f) -  bfImg.At(r,c,0);
        }
        mapping[count++] = r * raw->cols + c;
      }
    }
  }
  for (size_t r = 0; r < mTraces.n_rows; r++) {
    for (size_t c = 0; c < mTraces.n_cols; c++) {
      assert(isfinite(mTraces.at(r,c)));
    }
  }
  assert(mapping.size() == (size_t)count);

  /* Subtract off the median. */
  fmat colMed = median(mTraces);
  frowvec colMedV = colMed.row(0);
  for(size_t i = 0; i < mTraces.n_rows; i++) {
    mTraces.row(i) = mTraces.row(i) - colMedV;
  }

  /* Get the quantiles of the mean difference for well and exclude outliers */
  fmat mad = mean(mTraces, 1);
  fvec madV = mad.col(0);
  mMadSample.Clear();
  mMadSample.Init(1000);
  mMadSample.AddValues(madV.memptr(), madV.n_elem);
  float minVal = mMadSample.GetQuantile(.25) - iqrThreshold * mMadSample.GetIQR();
  float maxVal = mMadSample.GetQuantile(.75) + iqrThreshold * mMadSample.GetIQR();
  for (size_t i = 0; i < madV.n_rows; i++) {
    if (madV[i] <= minVal || madV[i] >= maxVal) {
      mWells[mapping[i]] = Filtered;
    }
  }
}