search

rbeyer / hiproc

Python libraries and programs for processing HiRISE data.
Apache License 2.0
4 stars 3 forks source link

Comparison of UA Dejitter HiRISE and HiPROC Dejitter HiRISE #6

Open jlaura opened 3 years ago

jlaura commented 3 years ago

Description

Seeing differences in the ephemeris information between UA dejittered, noproj data, and hiproc dejittered noproj data. Seeking input on the differences.

What I Did

Ran the full pipeline on ESP_018957_2655

using: hiproc --keep --verbose --precision --max_workers 8 "$@"

This ran HiJACK and generated a RED NOPROJ mosaic as expected (awesome). I then compared a UA dejittered, NOPROJ RED cube and the hiproc generated cube, looking at the ephemeris information.

Ephemeris times and look vectors are identical between the two data products. Again, awesome! The spacecraft XYZ positions differ though. isisA is the UA dejittered RED NOPROJ cube and isisB is the hiproc dejittered RED NOPROJ cube. Differences in spacecraft position are measures at 100 points over the image. I have truncated the results below. The interesting thing is that the differences in x, y, and z are relatively consistent along and across track. I say relatively because the magnitude of the difference is monotonically increasing across and along track.

Is this expected behavior? Where might this systematic offset be coming from? Operator error?

isisA pos x isisA pos y isisA pos z isisB pos x isisB pos y isisB pos z diffx diffy diffz magnitude
-286308.739429 81203.708950 3.682743e+06 -286295.489350 81209.926531 3.682744e+06 -13.250079 -6.217581 -0.860491 14.661628
-286242.472858 81226.509431 3.682748e+06 -286229.222798 81232.727130 3.682749e+06 -13.250060 -6.217699 -0.860207 14.661645
-286176.206308 81249.309623 3.682753e+06 -286162.956267 81255.527440 3.682754e+06 -13.250041 -6.217817 -0.859923 14.661661
-286109.939398 81272.109657 3.682758e+06 -286096.689376 81278.327592 3.682758e+06 -13.250022 -6.217935 -0.859639 14.661677
-286043.672318 81294.909469 3.682762e+06 -286030.422316 81301.127522 3.682763e+06 -13.250003 -6.218053 -0.859355 14.661693
... ... ... ... ... ... ... ... ... ...
-285977.405260 81317.708991 3.682767e+06 -285964.155277 81323.927162 3.682768e+06 -13.249984 -6.218171 -0.859072 14.661709
-285911.137842 81340.508356 3.682772e+06 -285897.887878 81346.726645 3.682772e+06 -13.249964 -6.218289 -0.858788 14.661726
-285844.870255 81363.307498 3.682776e+06 -285831.620310 81369.525905 3.682777e+06 -13.249945 -6.218408 -0.858504 14.661742
-285778.602690 81386.106351 3.682781e+06 -285765.352764 81392.324876 3.682782e+06 -13.249926 -6.218526 -0.858220 14.661758
-285712.334765 81408.905046 3.682786e+06 -285699.084857 81415.123690 3.682786e+06 -13.249907 -6.218644 -0.857936 14.661774

@rbeyer If you want the UA cub and the notebook I am using for comparisons just let me know an I can get you a URL / gist.

rbeyer commented 3 years ago

Hmmmmmmmmm. First off, I'm not 100% sure, but I have opinions (you know, as always).

I will say that we uncovered some bugs working together with the UofA in their resolve_jitter program in 2012. We worked on a conversion of their MATLAB version of resolve_jitter to C++ in order to free their processing nodes from needing MATLAB licenses and to speed up processing. At the time, they did not have the capacity to investigate these issues, critically evaluate our changes, or even change over to this code, and they continue to use the MATLAB version to this day, and incur the costs of licensing and extra processing time (which always made me sad).

Oleg's C++ version (which has not been publicly released) was much kinder and had #if statements that allowed you to compile the code with the bugs faithfully reproduced (bless him) or with the fixes enabled. I didn't bother with that nonsense in the Python version of resolve_jitter, I just fixed the damn bugs.

Anyway, resolve_jitter is what produces the files that are fed to hijitter to actually do the ephemeris corrections, and so some minor changes there could be the source of the differences you're seeing.

I suppose I could write the bugs into the Python version (please don't make me), or I could run hiproc using the "with bugs" version of Oleg's resolve_jitter on your image (may take some time to get to), and then you can compare the results maybe, and see if that's the difference?

I guess it is also possible that the SPICE data that were available when the U of A created this dejittered product have gotten a reconstruction since then, and this is entirely due to different SPICE kernels used? That one seems like a long shot. Not sure if those are carried along with the dejittered product you have or not so you could compare?

jlaura commented 3 years ago

I wonder if this is not, in some way, related that this ISIS issue then: https://github.com/USGS-Astrogeology/ISIS3/issues/3257. The offsets in that issue look like they are along track and of a lesser magnitude, but it would not be surprising if these are somehow related.

I also do not want to reproduce bugs in order to get a matching solution for the sake of matching! I realize too that I do not know when the image I tested against was dejittered, so I do not know what version of the UA pipeline was used. Glad to know a bit more of the lineage, but bummed to see another variable added. 😢

If you have some time this FY, it would be sweet to see if the 'with bugs' in Oleg's version produces output that is identical to the original UA output. If you are able to run it, feel free to just send me a link to the RED NOPROJ mosaic and I can download and run the comparison. No pressure, timeline, or rush getting to that as these offsets are quite small and quite consistent. I am going to proceed with testing these on ingest into Socet.

Thanks for the response and I'll keep you in the loop as I run these end-to-end tests (preempting your unit test style validation - sheesh, users are a pain! 😆).

rbeyer commented 3 years ago

I suppose it could be related to whatever the mystery is that is causing that issue. However, I hate to just wash my hands of it and say that "there was a change in the Matrix" and this is just where we are now, live with it. I will try and give this observation a run through using the version we have "with bugs" and let you know. Hopefully in the next month.

Of course, let's not forget that this could also be some wonderful new error that I have introduced by attempting to rewrite things.

jlaura commented 3 years ago

:+1: One more image to perhaps assist. Here is the diff image (using ISIS's fx). Left image is the diff, right image is the diff after running another round of processing that is done before ingest into socet.

image

The noise on the far left of the image is what I was hoping to see over the entire image (assuming the magnitude is low). The fact that we are seeing structure suggests to me that we have a systematic offset between the two results on some of the CCDs. Maybe? - I am not familiar enough with the sensor to do anything more than make a guess.

rbeyer commented 3 years ago

Hunh. So the Python HiCal program in hiproc contains improvements in noise processing that are currently staged for inclusion in the actual HiRISE pipelines. For this image (if I have it right in my head) that should only contribute an almost negligible downtrack gradient. There are more severe noise cases where new, extra processing kicks in, but it doesn't for this image.

What "structure" are you talking about? Are you talking about being able to see surface features in the left-hand diff, or are you talking about the vertical patterns (which are Channel and CCD joins) in the right-hand diff? And for whichever diff you're worried about, what are the stats? Whats the magnitude of the diff we're talking about here?

jlaura commented 2 years ago

@rbeyer Apologies for the delay getting back to you. The structure I am talking about are the surface features. I assumed that the vertical striping was a function of the CCD stitching. That vertical striping is visible in the left hand image two and I think that whatever processing we are doing to prepare for Socet is making that more pronounced in the right image.

Here are the stats for the left image:

Band 1 Block=500x500 Type=Float32, ColorInterp=Undefined
  Description = Red
  Minimum=-0.269, Maximum=0.154, Mean=-0.010, StdDev=0.002
  NoData Value=-3.40282265508890445e+38
  Mask Flags:
  Metadata:
    BANDWIDTH=300.000000
    BANDWIDTH_UNIT=NANOMETERS
    STATISTICS_MAXIMUM=0.15375693142414
    STATISTICS_MEAN=-0.0097139199355598
    STATISTICS_MINIMUM=-0.26859384775162
    STATISTICS_STDDEV=0.0023013041853695
    STATISTICS_VALID_PERCENT=99.09
    WAVELENGTH=700.000000
    WAVELENGTH_UNIT=NANOMETERS
rbeyer commented 2 years ago

Okay, so about the diff images: yeah, the vertical striping is from CCD and channel stitching artifacts, welcome to our HiRISE world of trying to make 20 different images "appear" to be just one image. I'm not sure what "additional" processing you're doing to get the right-hand diff, but whatever it is, it appears to be removing the surface feature structures, and while the vertical striping looks "more pronounced" in the right-hand diff it is probably the same, and just shows up better here, because the dynamic range has been collapsed by removal of the surface features, and qview auto-stretches (i.e., they were there, at that magnitude all along).

So I ran this observation through using the "with bugs" version of our C++ version of resolve_jitter and without, and then diff'ed the RED NOPROJ images, and the diff looks like your right-hand diff, no large image-scale surface features like you're seeing in the left-hand diff. If you zoom into the diff, you can see some small surface texture, but that's to be expected. So that means that the bug-fixes to resolve_jitter aren't the cause of the big differences you're seeing in the left-hand diff.

Unfortunately, this means that I don't know why the RED NOPROJ that hiproc is producing is so different from the one that you got from the UofA. Admittedly, I have not done a complete verification of all of the software pieces in hiproc against the output of the "living" HiRISE processing pipelines (I've just studied and aped them in Python code). So there is certainly the possibility that there are bugs that I have introduced in my implementation that could be the cause. When you see hiproc hit 1.0, then you'll know I've finished that verification. However, it is likely to take a long time, as even starting it is waiting on upstream changes in Tucson to settle down, and then I have to do a bunch of verification, which isn't as fun, so odds are good it'll get back-burnered.

Assuming that it isn't me (which, you know, it probably is), tracking down the source of the difference is going to be very difficult. We'd have to recapture the state of the HiRISE Team versions of their processing software at the time this cube you have was made, as well as the state of ISIS, and the ISISDATA area. It could be done, but it would be a lot of digital forensics work. That issue you referenced was about hundreds of HiRISE pixels of offset, and we still don't know why that is.

In an attempt to make us feel better about "letting go," the standard deviation of the diff is only 0.002 I/F, which tells us that the brightness variations between these two images are very, very small. Sure, there are outliers, but overall, maybe we can live with this?

Getting back to your original question: why are the reported spacecraft positions different by almost 15 m? That's got a more complicated answer that I don't know the answer to either (surprising, I know), but my guess is that it is related to the surface features we see between the images in your left-hand diff. The fact that we see structure that looks like surface features in the diff means that the surface features in the parent images aren't identically aligned in the two images, and since the NOPROJ images that you're seeing are the result of a dejittering process that is feature-based and outputs a new camera kernel, it is at least "consistent" such that seeing the structure differences in the left-hand diff should correlate with offsets in the spacecraft position. However, I don't know how to investigate that further, but I'm open to suggestions.

You should get an e-mail from NASA's Box service with a link to pull down the "with bugs" version of ESP_018957_2655_RED.NOPROJ.cub shortly.

jlaura commented 2 years ago

@rbeyer I spoke with the DTM group here and here is our plan given the challenges in 'ground' truthing this.

I am not 100% sure what the processing in the right image is. Something the DTM group does before ingestion into socet. I can inquire if you like? I am happy to simply let it go too. The post processing isn't something this library should ever worry about!

Thanks for processing through the C++ with bugs. It is a bummer to have the diffs, but at 0.002 I/F that is exceptionally small. The above proposed experiments are to test the livability of the approach. Before I say sure, I want the DTM generation folks to give the :+1: on the resultant DTM. So, tentative yes. I am really sensitive to the challenge here as well since 'truth' is not exactly a stable thing right now.

For the surface structure, I am pretty content to let that live for a bit to see what the resultant DTMs look like. We know that the HiRISE pipeline has reported similar shifts (albeit along track). It could be that the reported offsets are the egregious instances and that smaller offsets like these exist at different processing points. To really understand this, I suspect we would need a bunch of time and time machine.

So - very much appreciate you digging into this! Let me know what you think about the second item above, generating one more with bugs version from the C++ code. And we will go from there with generating a DTM and doing comparison of the results.

rbeyer commented 2 years ago

I am not 100% sure what the processing in the right image is. Something the DTM group does before ingestion into socet. I can inquire if you like? I am happy to simply let it go too. The post processing isn't something this library should ever worry about!

Nah, I don't care about the details, but the important thing is that after this processing the two versions of this "dejittered" cube are even more similar to one another, which means that "going into SOCET" there isn't much difference. So I will make the bold prediction that the resulting models will also be similar.

I am really sensitive to the challenge here as well since 'truth' is not exactly a stable thing right now.

For the surface structure, I am pretty content to let that live for a bit to see what the resultant DTMs look like. We know that the HiRISE pipeline has reported similar shifts (albeit along track). It could be that the reported offsets are the egregious instances and that smaller offsets like these exist at different processing points. To really understand this, I suspect we would need a bunch of time and time machine.

Yeah, this is a philosophy thing. We have a tendency to assume that something made by unknown other experts was "right" but in truth it could be just as flawed as something we make now.

And yes, short of a time machine, we may not even be able to per-pixel re-create the dejittered image that you have. So that might actually be another kind of test. Figure out how "old" your UofA dejittered cubes are, and if they're a few years old, then ask for a "new" version (maybe after HiRISE gets itself updated to ISIS 4.3 later this fall), and see if the new version is the same or different.

If it is different, then who knows what changed? If it is the same, then I've screwed up something somewhere in hiproc. I'll note that this kind of comparison is the sort of thing that I plan to do to get hiproc to version 1.0 as I mentioned above, I just haven't done it yet. So you could also wait on requesting this, as it kind of falls on my to-do list.

So - very much appreciate you digging into this! Let me know what you think about the second item above, generating one more with bugs version from the C++ code. And we will go from there with generating a DTM and doing comparison of the results.

I can do this, but I don't think there's much utility here, unless I'm misunderstanding something (which is always likely). The with-bugs-resolve_jitter version is very, very similar to the vanilla-hiproc version, and both of them are very different from the UofA version. Again, it is straightforward for me to do, but doing so implies a whole lot more work on your part downstream to make a third terrain model and do a bunch of comparison, and if we already know that the with-bugs-resolve_jitter version isn't substantially different from the vanilla-hiproc version, is there value in exercising that whole work process? Not sure. You tell me.

jlaura commented 2 years ago

I chatted some more with the DTM folks and we are good to go with the test using the dejittered data generated by this library and the older UA generated data. I'll keep you in the loop as they build the DTMs!