Camstats output doesn't make sense

[lwellerastro]

ISIS version(s) affected: 7.1.0, 7.2.0 I get expected/sensible output under 5.0.2, 6.0.0, 7.0.0

Description
Before recovering lost Kaguya TC data under 7.2.0, I compared the caminfo output of one image to that created under isis5.0.2 and found the camstats information to be quite different. cubediff run on the old image and this new image shows them to be identical. However, the caminfo/camstats reporting of emission and incidence angle information is quite different.

under 5.0.2 and 6.00     MinimumEmission          = 16.423659427042     MaximumEmission          = 19.293244403461     MinimumIncidence         = 59.443819396967     MaximumIncidence         = 61.109253930639

no caminfo output under 7.0.0 because of an introduced bug in caminfo

7.1.0 and 7.2.0 report     MinimumEmission                          = 1.0513035024623     MaximumEmission                          = 41.396795355947     MinimumIncidence                         = 35.210026880994     MaximumIncidence                         = 85.765737110116

These values make no sense for an equatorial image from a mapping mission. Additionally, basic geometry information in the image .lbl file includes INCIDENCE_ANGLE                      =  60.288 EMISSION_ANGLE                       =  16.580

which is similar 5.0.2 and 6.0.0 output. Note that the values for oblique resolution are also different between the versions and everything else is the same.

How to reproduce
See data under my work area Isis3Tests/Camstats/.
top level directory has output for isis7.2.0 Isis5.0.2/ contains data for that version for comparison Isis7.0.0 contains data for that version for comparison

Here's the processing history of the test image (img and lbl in directory for use). I removed path info for the shape model.

kaguyatc2isis from=TC1W2B0_01_00607S370E2634.lbl to=TC1W2B0_01_00607S370E2634.lev0.cub setnullrange=yes nullmin=-32768 nullmax=-0.00001 

spiceinit from=TC1W2B0_01_00607S370E2634.lev0.cub spksmithed=true shape=user model=Lunar_LRO_LOLAKaguya_DEMmerge_Global_512ppd_radius_demprep.cub

camstats from=TC1W2B0_01_00607S370E2634.lev0.cub attach=true linc=100 sinc=100 

footprintinit from=TC1W2B0_01_00607S370E2634.lev0.cub increaseprecision=FALSE inctype=vertices numvertices=40

caminfo from=TC1W2B0_01_00607S370E2634.lev0.cub to=TC1W2B0_01_00607S370E2634.lev0.pvl isislabel=true geometry=true originallabel=true statistics=true camstats=true linc=100 sinc=100 uselabel=true

I confirmed the differences exist at the camstats step when I output to a pvl file and make comparisons. Although caminfo/camstats was not generating output under 7.0.0, that version of camstats matches the output of isis5.0.2 which are the expected values.

I don't know if these differences are seen for other data sets.

Additional context

[Kelvinrr]

Can you post the kernels group for both images?

[lwellerastro]

isis7.2.0 image:

  Group = Kernels
    NaifCkCode                = -131350
    NaifFrameCode             = -131351
    LeapSecond                = $base/kernels/lsk/naif0012.tls
    TargetAttitudeShape       = ($base/kernels/pck/pck00009.tpc,
                                 $kaguya/kernels/pck/moon_080317.tf,
                                 $kaguya/kernels/pck/moon_assoc_me.tf)
    TargetPosition            = (Table,
                                 $kaguya/kernels/tspk/moon_pa_de421_1900-2050.-
                                 bpc, $kaguya/kernels/tspk/de421.bsp)
    InstrumentPointing        = (Table, $kaguya/kernels/ck/SEL_M_ALL_D_V02.BC,
                                 $kaguya/kernels/fk/SEL_V01.TF)
    Instrument                = $kaguya/kernels/ik/SEL_TC_V01.TI
    SpacecraftClock           = $kaguya/kernels/sclk/SEL_M_V01.TSC
    InstrumentPosition        = (Table,
                                 $kaguya/kernels/spk/SEL_M_071020_090610_SGMH_-
                                 02.BSP)
    InstrumentAddendum        = $kaguya/kernels/iak/kaguyaTcAddendum007.ti
    ShapeModel                = /work/projects/kaguya/TOPO/GlobalMerge/Lunar_-
                                LRO_LOLAKaguya_DEMmerge_Global_512ppd_radius_-
                                demprep.cub
    InstrumentPositionQuality = Reconstructed
    InstrumentPointingQuality = Reconstructed
    CameraVersion             = 2
    Source                    = isis
  End_Group

isis5.0.2 image:

  Group = Kernels
    NaifCkCode                = -131350
    NaifFrameCode             = -131351
    LeapSecond                = $base/kernels/lsk/naif0012.tls
    TargetAttitudeShape       = ($base/kernels/pck/pck00009.tpc,
                                 $kaguya/kernels/pck/moon_080317.tf,
                                 $kaguya/kernels/pck/moon_assoc_me.tf)
    TargetPosition            = (Table,
                                 $kaguya/kernels/tspk/moon_pa_de421_1900-2050.-
                                 bpc, $kaguya/kernels/tspk/de421.bsp)
    InstrumentPointing        = (Table, $kaguya/kernels/ck/SEL_M_ALL_D_V02.BC,
                                 $kaguya/kernels/fk/SEL_V01.TF)
    Instrument                = $kaguya/kernels/ik/SEL_TC_V01.TI
    SpacecraftClock           = $kaguya/kernels/sclk/SEL_M_V01.TSC
    InstrumentPosition        = (Table,
                                 $kaguya/kernels/spk/SEL_M_071020_090610_SGMH_-
                                 02.BSP)
    InstrumentAddendum        = $kaguya/kernels/iak/kaguyaTcAddendum007.ti
    ShapeModel                = /work/projects/kaguya/TOPO/GlobalMerge/Lunar_-
                                LRO_LOLAKaguya_DEMmerge_Global_512ppd_radius_-
                                demprep.cub
    InstrumentPositionQuality = Reconstructed
    InstrumentPointingQuality = Reconstructed
    CameraVersion             = 2
    Source                    = isis
  End_Group

I have also compared the output of catlab and cathist. The labels are identical and the only difference in history is the isis version.

[lwellerastro]

Some additional information - I have run phocube on the 5.0.2 and 7.2.0 cubes (includes the emission, incidence and phase angle backplanes) under their respective versions and the images are identical.

cubediff from=TC1W2B0_01_00607S370E2634.lev0.pho.cub from2=Isis5.0.2/TC1W2B0_01_00607S370E2634.lev0.pho.cub to=diff_phocube

  Group = Results
    Compare = Identical
  End_Group

The range of incidence and emission angles on the phocube file better matches the output reported via isis5.0.2 camstats when the program stats is run on file.

Incidence band stats output for phocube file: Average = 60.279095840589 StandardDeviation = 0.35080094295588 Variance = 0.12306130157874 Median = 60.280742770354 Mode = 60.429179303102 Skew = -0.014084310183994 Minimum = 59.443820953369 Maximum = 61.109252929688

Camstats Incidence angle information for isis5.0.2 image:

egrep Incidence 5.0.2_camerastats.pvl 
Group = IncidenceAngle
  IncidenceMinimum           = 59.443819396967
  IncidenceMaximum           = 61.109253930639
  IncidenceAverage           = 60.275539493619
  IncidenceStandardDeviation = 0.36112761738122

and isis7.2.0 image:

egrep Incidence 7.2.0_camerastats.pvl
Group = IncidenceAngle
  IncidenceMinimum           = 35.210026880994
  IncidenceMaximum           = 85.765737110116
  IncidenceAverage           = 59.29816852649
  IncidenceStandardDeviation = 9.1388846253812

I have found similar behaviour from camstats vs phocube for LROC NAC images while processing under 8.0.0-RC1.

[acpaquette]

I think this PR is likely the culprit for this change: https://github.com/DOI-USGS/ISIS3/pull/4600

With the change being suggested by Kris under this issue: https://github.com/DOI-USGS/ISIS3/issues/3600

[acpaquette]

@lwellerastro Can you run the same tests as above but with the shape set to an ellipsoid? The changes seen in the above linked PR likely make some significant differences if you are using a more detailed shape

[lwellerastro]

@acpaquette, was just doing that. The camstats output of those runs (spiceinit shape=ellipsoid) under 5.0.2 and 7.2.0 are identical except for the Oblique Resolution Standard Deviations and those diffs are in in the 13th decimal location. The values also make more sense for illumination angles for this type of data:

Group = EmissionAngle
  EmissionMinimum           = 16.453545888993
  EmissionMaximum           = 19.321301051227
  EmissionAverage           = 17.492384384575
  EmissionStandardDeviation = 0.90135523320707
End_Group

Group = IncidenceAngle
  IncidenceMinimum           = 59.440899879328
  IncidenceMaximum           = 61.139443826058
  IncidenceAverage           = 60.285995955023
  IncidenceStandardDeviation = 0.36934260759414
End_Group

I don't consider the DEM file I'm using to be especially detailed at 60 m/pixel for data that are generally around 10 m/pix. And this particular KaguyaTC image is near the equator (lat -37), covers about 35 km across and although not completely nadir, it's certainly not oblique. I wouldn't expect emission or incidence angles to vary that much across an image of this type as reported by the more current versions of camstats. I'd expect to see this sort of variation for flyby images having disk, limb and terminator in the scene. And I understand why irregular bodies might need local emission over the ellipsoid model, but I'm not sure it makes sense for mapping data.

The LROC NAC image (resolution 0.6 m/pixel) I happened to look at is very near the south pole (-86, but does not cross) and has high incidence angle for sure, but I find it hard to believe the range camstats reports (under 8.0.0RC1 using a 60m/pix dem) which is 69 to 114 degrees. That latter value seems way out of bounds for having an average emission angle of 13 (range is 0 to 36 which is very large for this sort of image; phocube says 0-3 which makes more sense)), essentially looking down. There is no "space" in the image, only heavy shadow and some illuminated features. I can't find that value via qviewing the image and using the tracking tool, but the illumination values constantly change from pixel to pixel so who knows. And tracking shows identical values for "incidence" vs "local incidence" if that means anything.

[acpaquette]

Thanks @lwellerastro I do think that provides some good context here. I think the best thing we could do is some manual verification of values by running campt and then computing the Emission, Incidence etc. angles manually and conferring with ISIS that those are the correct values

[acpaquette]

@lwellerastro Good news on this, the reason for this behavior has been found. Bad news is that we are going to have to wait until next support sprint to get the fix in. Basically after the addition of computing local photometric data there is a bug where the default emission angle computation uses local normal vector (dem normal vector) instead of the ellipsoid normal vector. This results in apps like camstats actually reporting the local emission angle rather than the ellipsoid emission angle.

If you would like to verify, run phocube with localemission and localincidence enabled, then run stats on the output cube. You'll find that the emission band looks fine, and the local emission band is producing nearly the same statistical values that camstats is producing.

Ellipsoid Incidence (band 2) and Emission (band 3):

Local (DEM) Incidence (band 4) and Emission (band 5):

In terms of a fix, the internal state between the DEM shape model and the ellipsoid model may need tweaking. I am going to get a bit into the weeds incase someone else picks this up other than myself.

The DEM and shapemodel share an internal variable that stores the normal vector (m_normal). If we compute the local photometric angles first, we store the local normal in this variable and set that a normal has been computed to true. If we then ask for the emission angle, incidence angle, etc. from the cube given a sample and line, those functions check if a normal has been computed. Since the local normal has been computed and stored in m_normal, we just end up using that local normal.

To some degree this may be considered desired behavior as Tim indicates in #4600, when he mentions that truth data for our tests will need to be updated. In that PR we did end up updating the camstats tests to account for this change.

It seems like some further discussion should be had here and that we need to get all of our apps inline so that they produce the same expected values. It seems like another possible solution would be for camstats to also potentially report local photometric angles as desired.

[lwellerastro]

Thanks for the detailed explanation @acpaquette. I understand what you are saying is the problem at the program level. But that very much changed the behavior and output of the program. I think an option for local angles needs to exist and it shouldn't be the default since it never was prior to the change.

I understand why some would say that intersecting a shape model and basing results on that is the desirable/only way to go, but I find the general camera statistics based on an ellipsoid indispensable for much of what I do. I am using the output of caminfo camera statistics in particular as input to database tables which when queried for emission and incidence help me determine if my data have a limb or terminator in the scene. The problem with intersecting the shape for that information is my equatorial, nadir images appear to be somewhere else because they happen to have strong topography that covers nearly the full spectrum of available emission and incidence values. I can't use the output anymore to make critical decisions about the data - for example, maybe some data are off nadir and having more meaningful emission angles would help identify those. I would prefer not having a second set of data that is initialized with spice based on an ellipsoid to work around the issue, especially when working with very large datasets.

I like the idea of basing the photometric angles (and radius) on the elliptical and producing the min/max values for the same based on local shape model feedback as an option for camstats in the same manner that phocube does.

[jlaura]

I would think that the return values would be dependent upon the shape model that was defined during spiceinit. If one initializes the cube with an ellipsoidal shape model, I would assume that photometric angles would be based on the ellipsoid, same if I initialize with a high resolution DTM. Then I would assume photometric angles are based on the high res DTM.

Is that what the current behavior is or is something different going on?

[lwellerastro]

@jlaura, the current behavior appears to be based on the shape model specified during spiceinit. This is a new behavior. The change seems to have occurred with isis7.1.0 based on testing. Prior to that, and for many years (decade at least), the photometric values were based upon an ellipsoid.

I immediately noticed the change when I had to reprocess Kaguya TC data due to a disk failure and processed under a more recent version of isis (data were originally processed under 5.0.2). The range of incidence angles was suddenly rather large for my data set and honestly didn't make sense for the images I was looking at (equatorial).

The only new things I've seen in regard to caminfo/camstats is being able to use the blob information already present on the image when generating that output for caminfo, and the post @acpaquette referenced in regard to using local emission for oblique data resolution logic in #3600 which unintentionally introduced the change.

[jlaura]

Sounds good. My preferred solution would be to have the returned values be based on the shape model being used. I get that this is new behavior. For anyone not accustomed to the old behavior this seems a lot more intuitive. The angles reported are based on the shape used. I would be all for adding an argument that let one override that default to say give me these angles but use the sphere/ellipsoid. This way, one wouldn't need to spiceinit twice depending on their use case.

I agree 100% that having the info on the ellipsoid is important for creating those database and building reasonable summary statistics. It is also important to get the right answer based on the shape used.

Glad that this came to light and I am looking forward to seeing how it is resolved.

Just my 2¢.

[lwellerastro]

If the current behavior is maintained and an option is added to get statistics based on an ellipsoid/sphere, then other programs may need to change as well to be consistent, something @acpaquette was pointing out in his last comments.

Phocube generates photometric backplanes based on an ellipsoid but has the option to also create local emission and incidence backplanes based on the shape model on the image label. Can't think of what else operates that way off the top of my head, but if a DEM would take precedence because the image is initialized with one (and this honestly does make a lot of sense), then other programs would need to be evaluated and possibly changed. Could be a bit of a project.

I suspect things operate the way they do simply because high quality shape models were not available for most bodies until recently. So it's not out of the question to discuss how this should be handled these days.

[KrisBecker]

Shape model interaction within the Camera system is complex. @lwellerastro is correct that it was introduced as a side effect of #3600, however it could also occur in other geometric calculations in the ISIS system. But how and why it is occurring is not apparent because of the complexity and interactions of Camera and ShapeModel classes. There are several solutions to the problem, but it really is an indication of a larger problem with the design of the ShapeModel class and all of its derivatives.

Essentially, the root cause of the problem is that the default behavior of the oblique detector resolution methods is to use the local normal. All other instances (via EmissionAngle()) use whatever normal currently exists if the derived shape model does not reimplement emissionAngle().

The explanation for why it impacts campt, caminfo and camstats is because the tools they each use, CameraPointInfo in campt, CamTools::BandGeometry() in caminfo, and CameraStatistics in camstats, all compute the oblique detector data before the emission angle, which means by the time the emission using EmissionAngle() and incidence (via ShapeModel::incidenceAngle()) angles are computed, it uses the preexisting local normal computed by the oblique detector methods. phocube appears to be self-aware of this issue.

Its not clear the best way to fix this problem, but some (arguably not so great) suggestions would be to move computation of oblique detector values after other geometric values are computed, change the default behavior of the oblique detector functions to use the ellipsoid, always force recalculation of the surface normal in all ShapeModel::emissionAngle(), or reimplement ShapeModel::emissionAngle() (and also ShapeModel::incidenceAngle()) in every class derived from ShapeModel (note EquatorialCylindricalShape and DemShape currently do not). All of these solutions have likely undesirable side effects and potential pitfalls.

It may be time to consider redesign of the ShapeModel class and potentially the Camera class.

[acpaquette]

That is the same conclusion I came to @KrisBecker minus the Camera class redesign. Also I am not sure what you mean when you say that phocub is self aware but it just does the computation in a different order. First computing the default photogrammetric angles then doing the local photogrammetric angles.

Having inconsistent state is a significant problem and we need to decide on some expected output. I think revisiting the ShapeModel stack of classes would be useful, potential removing some of the default behavior in the parent ShapeModel class and forcing the child classes to implement said behavior could be a fairly do-able solution in the short term. Anything beyond that will likely need to be it's own IAA which may not be possible at this point.

[acpaquette]

I also made a discussion post #5275 If we want to move over there to prevent overloading this issue with unnecessary context.