Current Situation

UnitMap starts with a dimensionless pixel and beam in the "Customary Map":

UnitMap::mapCust->insert(map<String, UnitName>::value_type("pixel", 
            UnitName("pixel", UnitVal(1.,UnitDim::Dnon), "pixel")));

and dimensionless beam and pixel can also be defined if FITSunits are loaded:
UnitName("BEAM", UnitVal(1.0, getStringFITS(0)), "dimensionless beam"), 
UnitName("PIXEL", UnitVal(1.0, getStringFITS(12)),"dimensionless pixel"),

The problem that provoked my inquiry:

In ComponentImager::project:

Unit fluxUnits;
{
Unit imageUnit = image.units();
const String& imageUnitName = imageUnit.getName();
UnitMap::putUser("pixel", UnitVal(pixelLatSize.radian() * pixelLongSize.radian(), "rad.rad"));
const Unit pixel("pixel");
if (imageUnitName.contains("pixel")) {
   // Get the new definition of the imageUnit which uses the new
   // definition of the pixels unit.
    imageUnit = image.units();
   fluxUnits.setValue(imageUnit.getValue() * pixel.getValue());
   fluxUnits.setName(imageUnitName + String(".") + pixel.getName()); }

[...]
const Unit jy("Jy");
if (fluxUnits != jy) { os << LogIO::WARN

What I think this is supposed to do, is
  1. define a "pixel" unit in the User Map, with the expectation that the User Map will be searched before the Custom one.
  2. creating the new Unit pixel("pixel") is supposed to clear the UnitCache and create a pixel unit with dimensions
  3. the second call to imageUnit = image.units() is supposed to be parsed using the updated UnitMap, i.e. in this case image.units()="Jy/pixel", so its supposed to turn that pixel into a solid angle.
  4. Then the solid angles cancel in fluxUnits, and all is well in the check == Jy.

But it doesn't work. the new Unit pixel does indeed have dimensions of solid angle, but imageUnit remains "1e-26 kg.s-2._-1" i.e. the dimension that's not a dimension "_" sticks around. No combination of deleting Jy/pixel from the Cache, moving the pixel from Custom to User Map, etc has succeeded in fixing this.

Generalization

The nondimensioned pixel can cause problems, because can end up in the Map and Cache multiple times, and different parts of the code use it differently, making it difficult to make any assumptions about how it will be defined at any given point.

An important example, and IMO an illustration of how pixel should be used is most of SkyComp and SkyCompRep e.g. in ::convertToJy and ::integralToPeakFlux: defineBrightnessUnits() defines the pixel and beam Units using the Coordsys and imageInfo, and when the method is done, undefineBrightnessUnits removes those from the UnitMap again. (note that this is one of the reasons you can't assume anything about the value of pixel right now)
  Unit brightnessUnits = defineBrightnessUnits(os, brightnessUnitOut, dirCoord,
                                               restoringBeam, integralInJy);
// Find peak value. Because we have defined /beam and /pixel units
// above we can use Quanta mathematics to get the answer we want

// /beam or /pixel units divided out here
   Quantum fluxPeak = fluxIntegral / tmp;    

   fluxPeak.convert(brightnessUnits);
   parameters(0) = fluxPeak.getValue();
   
// Undefine /beam and /pixel units
   SkyCompRep::undefineBrightnessUnits();

Proposed Solution 1: Don't use nondimensional "pixel" in the unitMap

Routines that load images would have to try to define pixel according to the Coordinatesystem, instead of making it nondimensional. I think if an image is trying to have units of flux/pixel, but doesn't have a valid coordinatesystem from which the pixel size in Sr can be determined, the brightness unit should stay undefined and a warning printed. Examples of current code that would need to be modified:
void HDF5Image::restoreUnits (const RecordInterface& rec)
void PagedImage::restoreUnits (const TableRecord& rec)

  if (! unitName.empty()) {
    // OK, non-empty unit, see if it's valid, if not try some known things to
    // make a valid unit out of it.
    if (! UnitVal::check(unitName)) {
      // Beam and Pixel are the most common undefined units
      UnitMap::putUser("Pixel",UnitVal(1.0),"Pixel unit");
      UnitMap::putUser("Beam",UnitVal(1.0),"Beam area");
    }

The /linear/ pixel unit is used in image analysis in ways such as this:
ImageProfileFit::ImageProfileFit(const ImageProfileFit& other)
   UnitMap::putUser("pix",UnitVal(1.0), "pixel units");
This is apparently so that the output of a fit can refer to an axis in "pix" units without having a coordsys.

Comments such as this one
   // This will probably be fixed in the future by doing the fits
   // in pixel space here and requiring the caller to deal with converting
   // to something astronomer friendly if it so desires.
Imply that this use of pixel and user units for spatial dimensions may be under revision.

In any case, I don't see a problem with using "pix" as a linear coordinate unit in this way, although it might be better to convert to something like "linearpix" to distinguish from "pixel" which is most commonly a solid angle.

String ImageFitter::_fluxToString(uInt compNumber)
may be making assumptions about pixel and beam, in particular it uses
 
       Quantity intensityToFluxConversion = peakIntensity.getUnit().contains("/beam")
                ? Quantity(1.0, "beam") : Quantity(1.0, "pixel");
so that Quantum math can do the peak to integrated flux conversions; I think that requiring that "pixel" have actual dimensions equivalent to Sr would not break this, and would in fact make the necessary conversions more robust, but it would need checking.

void Image2DConvolver::dealWithRestoringBeam
and Double ImageRegrid::findScaleFactor

deal extensively with Jy/pixel and Jy/beam but use Coordsys to do the pixel calculations, and not Quantum math directly, so I think would be ok.

What about Beam?

IMO one should never have an image whose units are /pixel without having a coordsys that defines pixel in Sr. If one does have such a thing, the brightness units are meaningless, since it'd be Jy per completely unknown area on the sky.

In most cases, images with units of Jy/beam have a defined beam, so the same effort should be made in loading those images to define beam properly in the UnitMap. However, I can see the utility of processing images in Jy/beam even if the beam is not (well) known. Some actions, such as the flux of an extended source, aren't possible, but many things are. I'm not sure what the best solution here is - probably the easiest is to allow a nondimensional beam in the UnitMap, which will cause the appropriate calculations to fail, and allow others to proceed.

In any case, there still needs to be a decision IMO about where beams and pixels should go - in the "custom" Map or the User Map (currently its both). If my interpretation of the code is correct, and something in the User Map is supposed to override the Cust Map, then usage such as in SkyCompRep where the User beam/pixel are defined and then undefined should be the recommended practice. However, if that overiding is supposed to happen its not clear to me why ComponentImager is broken smile

-- RemyIndebetouw - 2010-11-04
Topic revision: r3 - 2010-11-04, RemyIndebetouw
This site is powered by FoswikiCopyright © by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding NRAO Public Wiki? Send feedback