C++ KstVectorPtr::interpolate方法代码示例

//.........这里部分代码省略.........
    if (!yv->resize(_ns)) {
      // FIXME: handle error?
      unlockInputsAndOutputs();
      return false;
    }
    yv->zero();
    i0 = 0; // other vectors may have diffent lengths, so start over
    v_shift = _ns;
  } else {
    // calculate shift and new samples
    // only do shift optimization if all used vectors are same size and shift
    v_shift = (*_xInVector)->numShift();
    v_new = (*_xInVector)->numNew();

    for (KstVectorMap::ConstIterator i = VectorsUsed.begin(); i != VectorsUsed.end(); ++i) {
      if (v_shift != i.data()->numShift()) {
        v_shift = _ns;
      }
      if (v_new != i.data()->numNew()) {
        v_shift = _ns;
      }
      if (_ns != i.data()->length()) {
        v_shift = _ns;
      }
    }

    if (v_shift > _ns/2 || force) {
      i0 = 0;
      v_shift = _ns;
    } else {
      KstVectorPtr xv = *_xOutVector;
      KstVectorPtr yv = *_yOutVector;
      for (int i = v_shift; i < _ns; i++) {
        yv->value()[i - v_shift] = yv->value()[i];
        xv->value()[i - v_shift] = xv->value()[i];
      }
      i0 = _ns - v_shift;
    }
  }

  _numShifted = (*_yOutVector)->numShift() + v_shift;
  if (_numShifted > _ns) {
    _numShifted = _ns;
  }

  _numNew = _ns - i0 + (*_yOutVector)->numNew();
  if (_numNew > _ns) {
    _numNew = _ns;
  }

  (*_xOutVector)->setNewAndShift(_numNew, _numShifted);
  (*_yOutVector)->setNewAndShift(_numNew, _numShifted);

  double *rawxv = (*_xOutVector)->value();
  double *rawyv = (*_yOutVector)->value();
  KstVectorPtr iv = (*_xInVector);

  Equation::Context ctx;
  ctx.sampleCount = _ns;
  ctx.xVector = iv;

  if (!_pe) {
    if (_equation.isEmpty()) {
      unlockInputsAndOutputs();
      return true;
    }

    QMutexLocker ml(&Equation::mutex());
    yy_scan_string(_equation.latin1());
    int rc = yyparse();
    _pe = static_cast<Equation::Node*>(ParsedEquation);
    if (_pe && rc == 0) {
      Equation::FoldVisitor vis(&ctx, &_pe);
      KstStringMap sm;
      _pe->collectObjects(VectorsUsed, ScalarsUsed, sm);
      ParsedEquation = 0L;
    } else {
      delete (Equation::Node*)ParsedEquation;
      ParsedEquation = 0L;
      _pe = 0L;
      unlockInputsAndOutputs();
      return false;
    }
  }

  for (ctx.i = i0; ctx.i < _ns; ++ctx.i) {
    rawxv[ctx.i] = iv->value(ctx.i);
    ctx.x = iv->interpolate(ctx.i, _ns);
    rawyv[ctx.i] = _pe->value(&ctx);
  }

  if (!(*_xOutVector)->resize(iv->length())) {
    // FIXME: handle error?
    unlockInputsAndOutputs();
    return false;    
  }

  unlockInputsAndOutputs();
  return true;
}

void BinnedMap::binnedmap() {
  KstVectorPtr x = *_inputVectors.find(VECTOR_X);
  KstVectorPtr y = *_inputVectors.find(VECTOR_Y);
  KstVectorPtr z = *_inputVectors.find(VECTOR_Z);

  KstMatrixPtr map = *_outputMatrices.find(MAP);
  KstMatrixPtr hitsMap = *_outputMatrices.find(HITSMAP);

  if (autoBin()) {
    AutoSize(X(),Y(), &_nx, &_xMin, &_xMax, &_ny, &_yMin, &_yMax);
  }
  
  bool needsresize = false;
  if (_nx<2) {
    _nx = 2;
    needsresize = true;
  }
  if (_ny<2) {
    _ny = 2;
    needsresize = true;
  }
  
  if ((map->xNumSteps() != _nx) || (map->yNumSteps() != _ny) ||
       (map->minX() != _xMin) || (map->minY() != _yMin)) {
    needsresize = true;
  }
  
  if (map->xStepSize() != (_xMax - _xMin)/double(_nx-1)) {
    needsresize = true;
  }
  if (map->yStepSize() != (_yMax - _yMin)/double(_ny-1)) {
    needsresize = true;
  }
  
  if (needsresize) {
    map->change(map->tag(), _nx, _ny, _xMin, _yMin,
		(_xMax - _xMin)/double(_nx-1), (_yMax - _yMin)/double(_ny-1));
    map->resize(_nx, _ny);
    hitsMap->change(hitsMap->tag(), _nx, _ny, _xMin, _yMin,
		(_xMax - _xMin)/double(_nx-1), (_yMax - _yMin)/double(_ny-1));
    hitsMap->resize(_nx, _ny);
  }

  map->zero();
  hitsMap->zero();

  int ns = z->length(); // the z vector defines the number of points.
  double n,p, x0, y0, z0;
  for (int i=0; i<ns; i++) {
    x0 = x->interpolate(i, ns);
    y0 = y->interpolate(i, ns);
    z0 = z->interpolate(i, ns);
    p = map->value(x0, y0)+z0;
    map->setValue(x0, y0, p);
    n = hitsMap->value(x0, y0)+1;
    hitsMap->setValue(x0, y0, n);
  }

  for (int i=0; i<_nx; i++) {
    for (int j=0; j<_ny; j++) {
      p = map->valueRaw(i,j);
      n = hitsMap->valueRaw(i,j);
      if (n>0) {
	map->setValueRaw(i,j,p/n);
      } else {
	map->setValueRaw(i,j,KST::NOPOINT);
      }
    }
  }

  //calculate here...
}

KstObject::UpdateType KstPSDCurve::update(int update_counter) {
  int i_subset, i_samp;
  int n_subsets;
  int v_len;
  int copyLen;
  double mean;
  double y;
  bool force = false;
  KstVectorPtr iv = _inputVectors[INVECTOR];

  double *psd;

  if (KstObject::checkUpdateCounter(update_counter))
    return NO_CHANGE;

  if (update_counter <= 0) {
    force = true;
  } else {
    iv->update(update_counter);
  }

  v_len = iv->sampleCount();

  n_subsets = v_len/PSDLen+1;

  last_n_new += iv->numNew();

  if ((last_n_new < PSDLen/16) && (n_subsets - last_n_subsets < 1) && !force) {
    return NO_CHANGE;
  }

  psd = (*_sVector)->value();

  for (i_samp = 0; i_samp < PSDLen; i_samp++) {
    psd[i_samp] = 0;
  }

  for (i_subset = 0; i_subset < n_subsets; i_subset++) {
    /* copy each chunk into a[] and find mean */
    if (i_subset*PSDLen + ALen <= v_len) {
      copyLen = ALen;
    } else {
      copyLen = v_len - i_subset*PSDLen;
    }
    mean = 0;
    for (i_samp = 0; i_samp < copyLen; i_samp++) {
      mean += (
        a[i_samp] =
        iv->interpolate(i_samp + i_subset*PSDLen, v_len)
        );
    }
    if (copyLen>1) mean/=(double)copyLen;

    /* Remove Mean and apodize */
    if (removeMean() && appodize()) {
      for (i_samp=0; i_samp<copyLen; i_samp++) {
        a[i_samp]= (a[i_samp]-mean)*w[i_samp];
      }
    } else if (removeMean()) {
      for (i_samp=0; i_samp<copyLen; i_samp++) {
        a[i_samp] -= mean;
      }
    } else if (appodize()) {
      for (i_samp=0; i_samp<copyLen; i_samp++) {
        a[i_samp] *= w[i_samp];
      }
    }
    for (;i_samp < ALen; i_samp++) a[i_samp] = 0.0;

    /* fft a */
    rdft(ALen, 1, a);
    /* sum each bin into psd[] */
    psd[0]+=a[0];
    psd[PSDLen-1] += a[1];
    for (i_samp=1; i_samp<PSDLen-1; i_samp++) {
      psd[i_samp]+= cabs(a[i_samp*2], a[i_samp*2+1]);
    }
  }

  last_f0 = 0;
  last_n_subsets = n_subsets;
  last_n_new = 0;

  norm_factor = 1.0/(sqrt(double(Freq)*double(PSDLen))*double(n_subsets));

  psd[0]*=norm_factor;

  MaxY = MinY = mean = psd[0];
  if (psd[0]>0)
    MinPosY = psd[0];
  else (MinPosY = 1.0e300);

  /* normalize psd */
  for (i_samp=1; i_samp<PSDLen; i_samp++) {
    y = (psd[i_samp]*=norm_factor);
    if (y>MaxY)
      MaxY=y;
    if (y<MinY)
      MinY=y;
    if ((y>0) && (y<MinPosY))
//.........这里部分代码省略.........

void BinnedMap::binnedmap() {
    KstVectorPtr x = *_inputVectors.find(VECTOR_X);
    KstVectorPtr y = *_inputVectors.find(VECTOR_Y);
    KstVectorPtr z = *_inputVectors.find(VECTOR_Z);
    KstMatrixPtr map = *_outputMatrices.find(MAP);
    KstMatrixPtr hitsMap = *_outputMatrices.find(HITSMAP);
    KstScalarPtr autobin = *_inputScalars.find(AUTOBIN);

    if (autobin) {
        if (autobin->value() != 0.0) {
            _autoBin = true;
        } else {
            _autoBin = false;
        }
    }

    if (_autoBin) {
        double minx, miny, maxx, maxy;
        int nx, ny;

        autoSize(X(), Y(), &nx, &minx, &maxx, &ny, &miny, &maxy);

        setNX(nx);
        setNY(ny);
        setXMin(minx);
        setXMax(maxx);
        setYMin(miny);
        setYMax(maxy);
    } else {
        KstScalarPtr xmin = *_inputScalars.find(XMIN);
        KstScalarPtr xmax = *_inputScalars.find(XMAX);
        KstScalarPtr ymin = *_inputScalars.find(YMIN);
        KstScalarPtr ymax = *_inputScalars.find(YMAX);
        KstScalarPtr nx = *_inputScalars.find(NX);
        KstScalarPtr ny = *_inputScalars.find(NY);

        if (xmin) {
            _xMin = xmin->value();
        }
        if (xmax) {
            _xMax = xmax->value();
        }
        if (ymin) {
            _yMin = ymin->value();
        }
        if (ymax) {
            _yMax = ymax->value();
        }
        if (nx) {
            _nx = (int)nx->value();
        }
        if (ny) {
            _ny = (int)ny->value();
        }
    }

    bool needsresize = false;

    if (_nx < 2) {
        _nx = 2;
        needsresize = true;
    }
    if (_ny < 2) {
        _ny = 2;
        needsresize = true;
    }

    if ((map->xNumSteps() != _nx) || (map->yNumSteps() != _ny) ||
            (map->minX() != _xMin) || (map->minY() != _yMin)) {
        needsresize = true;
    }

    if (map->xStepSize() != (_xMax - _xMin)/double(_nx-1)) {
        needsresize = true;
    }
    if (map->yStepSize() != (_yMax - _yMin)/double(_ny-1)) {
        needsresize = true;
    }

    if (needsresize) {
        map->change(map->tag(), _nx, _ny, _xMin, _yMin,
                    (_xMax - _xMin)/double(_nx-1), (_yMax - _yMin)/double(_ny-1));
        map->resize(_nx, _ny);
        hitsMap->change(hitsMap->tag(), _nx, _ny, _xMin, _yMin,
                        (_xMax - _xMin)/double(_nx-1), (_yMax - _yMin)/double(_ny-1));
        hitsMap->resize(_nx, _ny);
    }

    map->zero();
    hitsMap->zero();

    int ns = z->length(); // the z vector defines the number of points.
    double n,p, x0, y0, z0;

    for (int i=0; i<ns; i++) {
        x0 = x->interpolate(i, ns);
        y0 = y->interpolate(i, ns);
        z0 = z->interpolate(i, ns);
        p = map->value(x0, y0)+z0;
        map->setValue(x0, y0, p);
//.........这里部分代码省略.........