Changeset 299


Ignore:
Timestamp:
01/26/05 02:00:06 (20 years ago)
Author:
kil064
Message:

add resample function

Location:
trunk/src
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/src/SDMath.cc

    r294 r299  
    599599  pTabOut->putSDHeader(sh);
    600600
    601 
    602601// Loop over rows and bin along channel axis
    603602 
     
    631630    pTabOut->putSDContainer(sc);
    632631  }
     632  return pTabOut;
     633}
     634
     635SDMemTable* SDMath::resample (const SDMemTable& in, const String& methodStr,
     636                              Float width) const
     637//
     638// Should add the possibility of width being specified in km/s. This means
     639// that for each freqID (SpectralCoordinate) we will need to convert to an
     640// average channel width (say at the reference pixel).  Then we would need 
     641// to be careful to make sure each spectrum (of different freqID)
     642// is the same length.
     643//
     644{
     645   Bool doVel = False;
     646
     647// Interpolation method
     648
     649  Int interpMethod = 0;
     650  convertInterpString(interpMethod, methodStr);
     651
     652// Make output table
     653
     654  SDMemTable* pTabOut = new SDMemTable(in, True);
     655
     656// Resample SpectralCoordinates (one per freqID)
     657
     658  const uInt nCoord = in.nCoordinates();
     659  Vector<Float> offset(1,0.0);
     660  Vector<Float> factors(1,1.0/width);
     661  Vector<Int> newShape;
     662  for (uInt j=0; j<in.nCoordinates(); ++j) {
     663    CoordinateSystem cSys;
     664    cSys.addCoordinate(in.getSpectralCoordinate(j));
     665    CoordinateSystem cSys2 = cSys.subImage(offset, factors, newShape);
     666    SpectralCoordinate sC = cSys2.spectralCoordinate(0);
     667//
     668    pTabOut->setCoordinate(sC, j);
     669  }
     670
     671// Get header
     672
     673  SDHeader sh = in.getSDHeader();
     674
     675// Generate resampling vectors
     676
     677  const uInt nChanIn = sh.nchan;
     678  Vector<Float> xIn(nChanIn);
     679  indgen(xIn);
     680//
     681  Int fac =  Int(nChanIn/width);
     682  Vector<Float> xOut(fac+10);          // 10 to be safe - resize later
     683  uInt i = 0;
     684  Float x = 0.0;
     685  Bool more = True;
     686  while (more) {
     687    xOut(i) = x;
     688//
     689    i++;
     690    x += width;
     691    if (x>nChanIn-1) more = False;
     692  }
     693  const uInt nChanOut = i;
     694  xOut.resize(nChanOut,True);
     695cerr << "width, shape in, out = " << width << ", " << nChanIn << ", " << nChanOut << endl;
     696//
     697  IPosition shapeIn(in.rowAsMaskedArray(0).shape());
     698  sh.nchan = nChanOut;
     699  pTabOut->putSDHeader(sh);
     700
     701// Loop over rows and resample along channel axis
     702
     703  Array<Float> valuesOut;
     704  Array<Bool> maskOut; 
     705  Array<Float> tSysOut;
     706  Array<Bool> tSysMaskIn(shapeIn,True);
     707  Array<Bool> tSysMaskOut;
     708  for (uInt i=0; i < in.nRow(); ++i) {
     709
     710// Get container
     711
     712     SDContainer sc = in.getSDContainer(i);
     713
     714// Get data and Tsys
     715   
     716     const Array<Float>& tSysIn = sc.getTsys();
     717     const MaskedArray<Float>& dataIn(in.rowAsMaskedArray(i));
     718     Array<Float> valuesIn = dataIn.getArray();
     719     Array<Bool> maskIn = dataIn.getMask();
     720
     721// Interpolate data
     722
     723     InterpolateArray1D<Float,Float>::interpolate(valuesOut, maskOut, xOut,
     724                                                  xIn, valuesIn, maskIn,
     725                                                  interpMethod, True, True);
     726     sc.resize(valuesOut.shape());
     727     putDataInSDC(sc, valuesOut, maskOut);
     728
     729// Interpolate TSys
     730
     731     InterpolateArray1D<Float,Float>::interpolate(tSysOut, tSysMaskOut, xOut,
     732                                                  xIn, tSysIn, tSysMaskIn,
     733                                                  interpMethod, True, True);
     734    sc.putTsys(tSysOut);
     735
     736// Put container in output
     737
     738    pTabOut->putSDContainer(sc);
     739  }
     740//
    633741  return pTabOut;
    634742}
     
    845953                           const casa::String& kernelType,
    846954                           casa::Float width, Bool doAll) const
     955//
     956// Should smooth TSys as well
     957//
    847958{
    848959
    849960// Number of channels
    850961
    851    const uInt chanAxis = asap::ChanAxis;  // Spectral axis
    852962   SDHeader sh = in.getSDHeader();
    853963   const uInt nChan = sh.nchan;
     
    8951005
    8961006    if (doAll) {
    897        uInt axis = asap::ChanAxis;
    898        VectorIterator<Float> itValues(valuesIn, axis);
    899        VectorIterator<Bool> itMask(maskIn, axis);
     1007       VectorIterator<Float> itValues(valuesIn, asap::ChanAxis);
     1008       VectorIterator<Bool> itMask(maskIn, asap::ChanAxis);
    9001009       while (!itValues.pastEnd()) {
    9011010
     
    9181027// Set multi-dim Vector shape
    9191028
    920        shapeOut(asap::ChanAxis) = valuesIn.shape()(chanAxis);
     1029       shapeOut(asap::ChanAxis) = valuesIn.shape()(asap::ChanAxis);
    9211030
    9221031// Stuff about with shapes so that we don't have conformance run-time errors
  • trunk/src/SDMath.h

    r294 r299  
    8888// Bin up spectra
    8989   SDMemTable* bin(const SDMemTable& in, casa::Int width) const;
     90
     91// Resample spectra
     92   SDMemTable* resample(const SDMemTable& in, const casa::String& method,
     93                        casa::Float factor) const;
    9094
    9195// Smooth
Note: See TracChangeset for help on using the changeset viewer.