NeoGeographyToolkit / StereoPipeline

The NASA Ames Stereo Pipeline is a suite of automated geodesy & stereogrammetry tools designed for processing planetary imagery captured from orbiting and landed robotic explorers on other planets.
Apache License 2.0
479 stars 169 forks source link

loading RPC model when stereo with Map-projected Images v2.6.1 #221

Closed yghlc closed 3 years ago

yghlc commented 5 years ago

Hi Thanks for providing such a useful tool. I'm processing images with RPC from Chinese satellites. The version is v2.6.1 with the binary packages on Github. I use this tool to orthorectify the images and extract DEM. But two issues confused me. (1) OSX and Linux have different behavior of copying RPC files. On OSX machine, mapproject would copy the RPC file without any change and set METADATATYPE=DG to the output file. On the Linux machine, this didn't happen. This is minor but makes me confusing on the difference between OSX and Linux version. I expected they were identical just running on different OS.

(2)the tool requires RPC files again when stereo with Map-projected Images. I followed the example of "Running Stereo with Map-projected Images". On Linux, the mapproject did not copy the RPC files, then the stereo tool complained RPCModel: GDAL resource appears not to have RPC metadata and crashed. I felt confused because the image already had been projected to an external DEM (WGS 84 / UTM zone 46N), why the tool required RPC file again? There is a similar question in Google Group: "correlation without loading camera model in v2.6.1" but haven't got a response.

Thanks for any suggestion or help.

Best Regards, Lingcao

oleg-alexandrov commented 5 years ago

Thank you for your feedback.

We currently do not have a Mac machine to verify what you say due to the government shutdown. We will, at some point. The tools are meant to work the same on Linux and Mac.

I will suggest you try using the daily build, at https://byss.arc.nasa.gov/stereopipeline/daily_build/ for Linux. I recall we fixed something regarding RPC recently.

The tools need the the RPC file again because when we do stereo, at triangulation phase, we need to go back to the original cameras, and those are not stored in the mapprojected image.

I did not respond yet to the other user's question because here it is weekend, and we are also under the government shutdown and we can't work. :) I will soon.

There is a good chance something is slightly different for your Chinese satellite which we never ran into. I will suggest you look at our manual more, particularly where we process GeoEye and CartoSat. Maybe what is there will help you.

You can also try to first run stereo without mapprojected images, just with the original images.

If you are still out of luck, and if you can share the images and cameras, I can take a look.

yghlc commented 5 years ago

@oleg-alexandrov Thank you very much for your response. I manually copied the RPC file, then the tool worked and outputted a reasonable DEM. But I haven't got time to validate the quality of the DEM. I will try your suggestions in the next few days and let you know the results. I may also look into the source codes if I have time. Sorry for that I cannot share the original data to you because of the agreement with the vendor.

oleg-alexandrov commented 5 years ago

I will also suggest you try parallel_stereo with --stereo-algorithm 2. That implements MGM, a newer algorithm, which can give much nicer results.

yghlc commented 5 years ago

@oleg-alexandrov I have tried your suggestions.

  1. I tried to run stereo without mapprojected images, with "-t rpc" session. Stereo exited with errors: *.tif contains a non-normal georeference. --> Inputs are map projected --> Generating a 3D point cloud. terminate called after throwing an instance of 'vw::cartography::ProjectionErr' what(): Bad projection in GeoReference.cc. Proj.4 error: latitude or longitude exceeded limits. This could due to the RPC file of the Chinese satellite. I opened the images using ENVI (both ENVI4.8 and ENVI 5.3). Only in ENVI 4.8, "Open External File->IKONOS->Open with RPC positioning" shown the correct coordinates. Some other information is in an XML but ASP cannot parse it. Interestingly, mapproject can correctly utilize the rpc file and output an orthorectified image using an external DEM. I downloaded the ASP codes but need more time to explore.

  2. The latest built ASP (Jan 19) on Linux don't copy the rpc files when running mapproject. So I wrote a script to copy them.

  3. I tried parallel_stereo with --stereo-algorithm 2 for the images. It took a longer time. And the result is wrong. See the figure below, image

Thanks again for the useful tool. I have a emerengy project need to work on recently. Hopefully, I will come back someday.

oleg-alexandrov commented 5 years ago

I appreciate you taking the time to play with our tools. Often the quality depends on the inputs, both the images and the cameras. Given that we never tried your satellite, I can't say more at this moment. In case you are able to specify the satellite name, and perhaps point us to some sample public images and cameras, in the future we can try to see what is going on and we could add support for it.

yghlc commented 5 years ago

Oh, sorry for forgetting to mention the satellite name, it is zhiyuan-3 (ZY-3). From I know, it seems that the original data are not publicly available. I will send you the data if I can find some sample public images and cameras. Thanks!

oleg-alexandrov commented 5 years ago

I was able to fetch a ZY-3 image and reproduce the problem. It was because we were expecting a file ending in _RPC.TXT but that satellite was providing _rpc.txt.

oleg-alexandrov commented 5 years ago

The fix should be in at https://byss.arc.nasa.gov/stereopipeline/daily_build/ in a day or two.

zzqstar commented 3 years ago

Sorry to distrub you, I have a similar question about using stereo command in Ames Stereo Pipline (ASP), that is, the "*_PRC.txt" is not supported by the command the "stereo" or "parallel_stereo" in ASP, how can I use the satellite images with txt-rpc files to generate the digital surface model? I am trying to writing the rpc files into the header information of satellite image. If it doesn't work, I will convert the the formats of rpc from txt into xml.

oleg-alexandrov commented 3 years ago

Above you wrote PRC.txt. I think it should be RPC.txt. Maybe that's the problem. This should work.

In either case, we use GDAL to load RPC, which finds the _rpc.txt or _RPC.txt, and this functionality has been there for a while, I doubt it is not working.

oleg-alexandrov commented 3 years ago

I think you should run stereo with just the image files, but the directory having those files must have _RPC.txt files which are loaded quietly.

zzqstar commented 3 years ago

Above you wrote PRC.txt. I think it should be RPC.txt. Maybe that's the problem. This should work.

In either case, we use GDAL to load RPC, which finds the _rpc.txt or _RPC.txt, and this functionality has been there for a while, I doubt it is not working.

Thank you very much, good work you have done! The problem was solved.

zzqstar commented 3 years ago

Sorry to distrub you again. I found a problem while using the ame stereo pipline. Two different terrain digital surface model were generated using the same data and same commands (again). The number of the point cloud varied while excuting the command again. There was also a slight overall shift between the two generated terrain models. You may have a check. The realted files and commands are listed here. The data download link is [https://www.isprs.org/data/ikonos/StereoMP_1m_BW_8bit.zip] The stereo.default and the commands are as follows, stereo.default: stereo-algorithm 1 cost-mode 4 corr-kernel 7 7 subpixel-mode 12 median-filter-size 3 alignment-method affineepipolar

Commands:

!/bin/bash

outputStereoPrefix='run_stereo/run' tgtGSD=0.6 searchradiusfactor=3 BWDImage='po_97258_pan_0000000.tif' FWDImage='po_97258_pan_0010000.tif' bundle_adjust './'$BWDImage './'$FWDImage -o run_bundle/run stereo_def='stereo.default' parallel_stereo -t rpc './'$BWDImage './'$FWDImage -s './'$stereo_def './'$outputStereoPrefix --bundle-adjust-prefix \ run_bundle/run point2dem --search-radius-factor $searchradiusfactor --tr $tgtGSD './'$outputStereoPrefix-PC.tif point2las --compressed -r Earth './'$outputStereoPrefix-PC.tif

oleg-alexandrov commented 3 years ago

I don't recall now if --stereo-algorithm 1 (SGM) is deterministic. I think --stereo-algorithm 0 (regular block matching) is.

A little variation is maybe expected, perhaps there are some random numbers somewhere. I don't think there should be a systematic shift.

I will also suggest you try --stereo-algorithm 2, as that one produces more pleasing output as far as I know.

ScottMcMichael commented 3 years ago

Bundle adjustment at least has a random element during IP matching so that could be the cause of the different results you are seeing. Try looking through the output logs from the two runs and try to find the point where they diverge. You can also try map projecting portions of the input images after bundle adjustment, if they don't line up then the output DEMs probably won't either. The download link you posted is not working for me so I was unable to check on it.

Scott


From: Oleg Alexandrov notifications@github.com Sent: Monday, January 18, 2021 10:08 AM To: NeoGeographyToolkit/StereoPipeline StereoPipeline@noreply.github.com Cc: Subscribed subscribed@noreply.github.com Subject: [External] Re: [NeoGeographyToolkit/StereoPipeline] loading RPC model when stereo with Map-projected Images v2.6.1 (#221)

I don't recall now if --stereo-algorithm 1 (SGM) is deterministic. I think --stereo-algorithm 0 (regular block matching) is.

A little variation is maybe expected, perhaps there are some random numbers somewhere. I don't think there should be a systematic shift.

I will also suggest you try --stereo-algorithm 2, as that one produces more pleasing output as far as I know.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHubhttps://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fprotect2.fireeye.com%2Fv1%2Furl%3Fk%3D9d26377f-c2bd0f89-9d2646d3-8627ced94d24-23acffdc40f9b5b4%26q%3D1%26e%3Df075c292-9f2d-4db1-8736-e7625f8f838c%26u%3Dhttps%253A%252F%252Fgithub.com%252FNeoGeographyToolkit%252FStereoPipeline%252Fissues%252F221%2523issuecomment-762401072&data=04%7C01%7Cscott.t.mcmichael%40nasa.gov%7C5c3433b9bb994f7a662e08d8bbdc14aa%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C637465901204081045%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=4hqFiy4psIwkGBFL4L83JrX2rnyyULYFAk24RUxDZAc%3D&reserved=0, or unsubscribehttps://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fprotect2.fireeye.com%2Fv1%2Furl%3Fk%3D63ce315a-3c5509ac-63ce40f6-8627ced94d24-e28e58c9dd7c7eb2%26q%3D1%26e%3Df075c292-9f2d-4db1-8736-e7625f8f838c%26u%3Dhttps%253A%252F%252Fgithub.com%252Fnotifications%252Funsubscribe-auth%252FABDVLRAKBK75ESLEUKQDLL3S2R2ITANCNFSM4GRHMLEA&data=04%7C01%7Cscott.t.mcmichael%40nasa.gov%7C5c3433b9bb994f7a662e08d8bbdc14aa%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C637465901204090981%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=SbkpVVEB2VyChmiw9jgUO8RK5lnnzS7tr2rz%2Bgh0R9U%3D&reserved=0.

This e-mail, including any attached files, may contain confidential information, privileged information and/or trade secrets for the sole use of the intended recipient. Any review, use, distribution, or disclosure by others is strictly prohibited. If you are not the intended recipient (or authorized to receive information for the intended recipient), please contact the sender by reply e-mail and delete all copies of this message.

zzqstar commented 3 years ago

Thank you for the reply. I checked again on the results and found that:

  1. The --stereo-algorithm 1 (SGM) is deterministic. Based on the same bundle adjustment result, two terrain moedels generated by --stereo-algorithm 1 (SGM) were exact the same;
  2. The results of the two bundle adjustments were not identical with same configuration and commands. This has contributed to the variation of the terrain models generated by Ames Stereo Pipline.
  3. The terrain models (point clouds in 2 ) got a change in Y direction.
  4. The data download link is https://www.isprs.org/data/ikonos . In this page, the data is under the "Stereo Map Projection Product". You may also Google with "isprs ikonos data set", and the data can be downloaded from the first item in the page.
  5. A small suggestion. You may have bundle adjustment updatd by eleminated the the random element during IP matching. The random element may cause a unstable stereo matching result. The random element can be replaced by a calculated result from the stereo pair.

Thank you. Zhiqiang zuo

adehecq commented 3 years ago

Hi, Yes, as you found out, bundle_adjust can sometime give different results. In particular, the initial IP matching can contain many outliers that need to be properly filtered out and will impact your output camera model. Can you maybe send the logs of your bundle_adjust calls? You should first check that several passes are run (--num-passes, the default should be 2 passes). This should normally get rid of the outliers. You could try to increase this to 3 or slightly change the IP filtering options like --remove-outliers-params, --epipolar-threshold, --ip-inlier-factor, --ip-uniqueness-threshold etc. Each time, you can check the output match files using stereo_gui to see if they make sense. Maybe you can find some settings that give you more consistent results. You can find all options of bundle_adjust here: https://stereopipeline.readthedocs.io/en/latest/tools/bundle_adjust.html

Amaury

zzqstar commented 3 years ago

Thank you for the reply. I have tried again the bundle adjustment and the parameters were changed: " bundle_adjust './'$BWDImage './'$FWDImage --num-passes 3 --ip-num-ransac-iterations 1000 --remove-outliers-params '75 2.5 1.5 2.5' -o run_bundle/run " The other commands were remain unchaged. There is still a shift between the two terrain models generated from runing the ASP twice. The logs of bundle_adjustment are as follows:

Logs of the bundle_adjustment at first time

/home/zzqstar/Desktop/StereoPipeline270/libexec/bundle_adjust ./po_97258_pan_0000000.tif ./po_97258_pan_0010000.tif --num-passes 3 --ip-num-ransac-iterations 1000 --remove-outliers-params 75 2.5 1.5 2.5 -o run_bundle/run

Vision Workbench log started at 2021-01-19 17:53:56.

2021-01-19 17:53:56 {0} [ console ] : Using session: rpc. 2021-01-19 17:53:56 {0} [ console ] : --> Computing statistics for ./po_97258_pan_0000000.tif 2021-01-19 17:53:56 {0} [ console ] : using downsample scale: 3 2021-01-19 17:53:56 {0} [ console ] : Writing stats file: run_bundle/run-po_97258_pan_0000000-stats.tif 2021-01-19 17:53:56 {0} [ console ] : ./po_97258_pan_0000000.tif: [ lo: 0 hi: 255 mean: 73.9326 std_dev: 54.3488 ] 2021-01-19 17:53:56 {0} [ console ] : Using session: rpc. 2021-01-19 17:53:56 {0} [ console ] : --> Computing statistics for ./po_97258_pan_0010000.tif 2021-01-19 17:53:56 {0} [ console ] : using downsample scale: 3 2021-01-19 17:53:56 {0} [ console ] : Writing stats file: run_bundle/run-po_97258_pan_0010000-stats.tif 2021-01-19 17:53:56 {0} [ console ] : ./po_97258_pan_0010000.tif: [ lo: 0 hi: 255 mean: 72.0923 std_dev: 52.0436 ] 2021-01-19 17:53:56 {0} [ console ] : Using session: rpc. 2021-01-19 17:53:56 {0} [ console ] : Loading camera model: ./po_97258_pan_0000000.tif 2021-01-19 17:53:56 {0} [ console ] : Using session: rpc. 2021-01-19 17:53:56 {0} [ console ] : Loading camera model: ./po_97258_pan_0010000.tif 2021-01-19 17:53:56 {0} [ console ] : Using session: rpc. 2021-01-19 17:53:56 {0} [ console ] : Using session: rpc. 2021-01-19 17:53:56 {0} [ console ] : --> Computing statistics for ./po_97258_pan_0000000.tif 2021-01-19 17:53:56 {0} [ console ] : --> Reading statistics from file run_bundle/run-po_97258_pan_0000000-stats.tif 2021-01-19 17:53:56 {0} [ console ] : ./po_97258_pan_0000000.tif: [ lo: 0 hi: 255 mean: 73.9326 std_dev: 54.3488 ] 2021-01-19 17:53:56 {0} [ console ] : --> Computing statistics for ./po_97258_pan_0010000.tif 2021-01-19 17:53:56 {0} [ console ] : --> Reading statistics from file run_bundle/run-po_97258_pan_0010000-stats.tif 2021-01-19 17:53:56 {0} [ console ] : ./po_97258_pan_0010000.tif: [ lo: 0 hi: 255 mean: 72.0923 std_dev: 52.0436 ] 2021-01-19 17:53:56 {0} [ console ] : --> Matching interest points in StereoSession. 2021-01-19 17:53:56 {0} [ console ] : Using epipolar threshold = 188.656 2021-01-19 17:53:56 {0} [ console ] : IP uniqueness threshold = 0.8 2021-01-19 17:53:56 {0} [ console ] : Datum: Geodetic Datum --> Name: WGS_1984 Spheroid: WGS 84 Semi-major axis: 6378137 Semi-minor axis: 6356752.3142451793 Meridian: Greenwich at 0 Proj4 Str: +proj=longlat +datum=WGS84 +no_defs 2021-01-19 17:53:56 {0} [ console ] : Skipping rough homography. 2021-01-19 17:53:56 {0} [ console ] : Looking for IP in left image... 2021-01-19 17:53:56 {0} [ console ] : Using 1309 interest points per tile (1024^2 px). 2021-01-19 17:53:56 {0} [ console ] : Detecting IP 2021-01-19 17:53:56 {0} [ console ] : Building descriptors 2021-01-19 17:53:56 {0} [ console ] : Found interest points: 4999 2021-01-19 17:53:56 {0} [ console ] : Recording interest points to file: run_bundle/run-po_97258_pan_0000000.vwip 2021-01-19 17:53:56 {0} [ console ] : Looking for IP in right image... 2021-01-19 17:53:56 {0} [ console ] : Using 1309 interest points per tile (1024^2 px). 2021-01-19 17:53:56 {0} [ console ] : Detecting IP 2021-01-19 17:53:57 {0} [ console ] : Building descriptors 2021-01-19 17:53:57 {0} [ console ] : Found interest points: 4999 2021-01-19 17:53:57 {0} [ console ] : Recording interest points to file: run_bundle/run-po_97258_pan_0010000.vwip 2021-01-19 17:53:57 {0} [ console ] : --> Matching interest points using the epipolar line. 2021-01-19 17:53:57 {0} [ console ] : Uniqueness threshold: 0.8 2021-01-19 17:53:57 {0} [ console ] : Epipolar threshold: 188.656 2021-01-19 17:53:57 {0} [ console ] : Matching forward 2021-01-19 17:53:58 {0} [ console ] : ---> Obtained 4999 matches. 2021-01-19 17:53:58 {0} [ console ] : Matching backward 2021-01-19 17:53:58 {0} [ console ] : ---> Obtained 4999 matches. 2021-01-19 17:53:58 {0} [ console ] : Matched 1085 points. 2021-01-19 17:53:59 {0} [ console ] : Removed 31 points in stddev filtering. 2021-01-19 17:53:59 {0} [ console ] : Reduced matches to 1054 2021-01-19 17:53:59 {0} [ console ] : * Writing match file: run_bundle/run-po_97258_pan_0000000po_97258_pan_0010000.match 2021-01-19 17:53:59 {0} [ console ] : IP coverage fraction = 1 2021-01-19 17:53:59 {0} [ console ] : Number of matches in run_bundle/run-po_97258_pan_0000000__po_97258_pan_0010000.match 1054 2021-01-19 17:53:59 {0} [ console ] : Building the control network for run_bundle/run-po_97258_pan_0000000po_97258_pan_0010000.match. 2021-01-19 17:53:59 {0} [ console ] : Loading GCP files... 2021-01-19 17:53:59 {0} [ console ] : --> Bundle adjust pass: 0 2021-01-19 17:53:59 {0} [ console ] : Writing initial condition files... 2021-01-19 17:53:59 {0} [ console ] : Starting the Ceres optimizer... 2021-01-19 17:54:00 {0} [ console ] : Solver Summary (v 1.14.0-eigen-(3.3.7)-lapack-suitesparse-(5.6.0)-cxsparse-(3.2.0)-eigensparse-openmp-no_tbb)

                                 Original                  Reduced

Parameter blocks 1055 1055 Parameters 3171 3171 Residual blocks 2108 2108 Residuals 4224 4224

Minimizer TRUST_REGION

Dense linear algebra library EIGEN Trust region strategy LEVENBERG_MARQUARDT

                                    Given                     Used

Linear solver DENSE_SCHUR DENSE_SCHUR Threads 4 4 Linear solver ordering AUTOMATIC 1053,2 Schur structure 2,3,6 2,3,6

Cost: Initial 8.403343e+01 Final 2.001809e+01 Change 6.401534e+01

Minimizer iterations 25 Successful steps 19 Unsuccessful steps 6

Time (in seconds): Preprocessor 0.001707

Residual only evaluation 0.072870 (25) Jacobian & residual evaluation 0.951939 (19) Linear solver 0.053973 (25) Minimizer 1.087919

Postprocessor 0.000118 Total 1.089744

Termination: CONVERGENCE (Parameter tolerance reached. Relative step_norm: 9.589535e-09 <= 1.000000e-08.)

2021-01-19 17:54:00 {0} [ console ] : Writing final condition log files... 2021-01-19 17:54:00 {0} [ console ] : Removing pixel outliers in preparation for another solver attempt. 2021-01-19 17:54:00 {0} [ console ] : Outlier statistics: b = -0.132968, e = 0.247976. 2021-01-19 17:54:00 {0} [ console ] : Removing as outliers points with mean reprojection error > 1.5. 2021-01-19 17:54:00 {0} [ console ] : Removed 1 outliers by reprojection error. 2021-01-19 17:54:00 {0} [ console ] : Removed 0 outliers based on disparity of ip. 2021-01-19 17:54:00 {0} [ console ] : IP coverage fraction after cleaning = 1 2021-01-19 17:54:00 {0} [ console ] : Saving 1052 filtered interest points. 2021-01-19 17:54:00 {0} [ console ] : Writing: run_bundle/run-po_97258_pan_0000000__po_97258_pan_0010000-clean.match 2021-01-19 17:54:00 {0} [ console ] : --> Bundle adjust pass: 1 2021-01-19 17:54:00 {0} [ console ] : Starting the Ceres optimizer... 2021-01-19 17:54:01 {0} [ console ] : Solver Summary (v 1.14.0-eigen-(3.3.7)-lapack-suitesparse-(5.6.0)-cxsparse-(3.2.0)-eigensparse-openmp-no_tbb)

                                 Original                  Reduced

Parameter blocks 1054 1054 Parameters 3168 3168 Residual blocks 2106 2106 Residuals 4220 4220

Minimizer TRUST_REGION

Dense linear algebra library EIGEN Trust region strategy LEVENBERG_MARQUARDT

                                    Given                     Used

Linear solver DENSE_SCHUR DENSE_SCHUR Threads 4 4 Linear solver ordering AUTOMATIC 1052,2 Schur structure 2,3,6 2,3,6

Cost: Initial 8.324747e+01 Final 1.950171e+01 Change 6.374576e+01

Minimizer iterations 16 Successful steps 11 Unsuccessful steps 5

Time (in seconds): Preprocessor 0.001503

Residual only evaluation 0.043898 (16) Jacobian & residual evaluation 0.477924 (11) Linear solver 0.027145 (16) Minimizer 0.553203

Postprocessor 0.000073 Total 0.554779

Termination: CONVERGENCE (Parameter tolerance reached. Relative step_norm: 2.923126e-09 <= 1.000000e-08.)

2021-01-19 17:54:01 {0} [ console ] : Writing final condition log files... 2021-01-19 17:54:01 {0} [ console ] : Removing pixel outliers in preparation for another solver attempt. 2021-01-19 17:54:01 {0} [ console ] : Outlier statistics: b = -0.133039, e = 0.247879. 2021-01-19 17:54:01 {0} [ console ] : Removing as outliers points with mean reprojection error > 1.5. 2021-01-19 17:54:01 {0} [ console ] : Removed 0 outliers by reprojection error. 2021-01-19 17:54:01 {0} [ console ] : Removed 0 outliers based on disparity of ip. 2021-01-19 17:54:01 {0} [ console ] : IP coverage fraction after cleaning = 1 2021-01-19 17:54:01 {0} [ console ] : Saving 1052 filtered interest points. 2021-01-19 17:54:01 {0} [ console ] : Writing: run_bundle/run-po_97258_pan_0000000__po_97258_pan_0010000-clean.match 2021-01-19 17:54:01 {0} [ console ] : No new outliers removed, and the algorithm converged. No more passes are needed. 2021-01-19 17:54:01 {0} [ console ] : Writing: run_bundle/run-po_97258_pan_0000000.adjust 2021-01-19 17:54:01 {0} [ console ] : Writing: run_bundle/run-po_97258_pan_0010000.adjust

end

Logs of the bundle_adjustment at second time

/home/zzqstar/Desktop/StereoPipeline270/libexec/bundle_adjust ./po_97258_pan_0000000.tif ./po_97258_pan_0010000.tif --num-passes 3 --ip-num-ransac-iterations 1000 --remove-outliers-params 75 2.5 1.5 2.5 -o run_bundle/run

Vision Workbench log started at 2021-01-19 17:54:05.

2021-01-19 17:54:05 {0} [ console ] : Using session: rpc. 2021-01-19 17:54:05 {0} [ console ] : --> Computing statistics for ./po_97258_pan_0000000.tif 2021-01-19 17:54:05 {0} [ console ] : using downsample scale: 3 2021-01-19 17:54:05 {0} [ console ] : Writing stats file: run_bundle/run-po_97258_pan_0000000-stats.tif 2021-01-19 17:54:05 {0} [ console ] : ./po_97258_pan_0000000.tif: [ lo: 0 hi: 255 mean: 73.9326 std_dev: 54.3488 ] 2021-01-19 17:54:05 {0} [ console ] : Using session: rpc. 2021-01-19 17:54:05 {0} [ console ] : --> Computing statistics for ./po_97258_pan_0010000.tif 2021-01-19 17:54:05 {0} [ console ] : using downsample scale: 3 2021-01-19 17:54:05 {0} [ console ] : Writing stats file: run_bundle/run-po_97258_pan_0010000-stats.tif 2021-01-19 17:54:05 {0} [ console ] : ./po_97258_pan_0010000.tif: [ lo: 0 hi: 255 mean: 72.0923 std_dev: 52.0436 ] 2021-01-19 17:54:05 {0} [ console ] : Using session: rpc. 2021-01-19 17:54:05 {0} [ console ] : Loading camera model: ./po_97258_pan_0000000.tif 2021-01-19 17:54:05 {0} [ console ] : Using session: rpc. 2021-01-19 17:54:05 {0} [ console ] : Loading camera model: ./po_97258_pan_0010000.tif 2021-01-19 17:54:05 {0} [ console ] : Using session: rpc. 2021-01-19 17:54:05 {0} [ console ] : Using session: rpc. 2021-01-19 17:54:05 {0} [ console ] : --> Computing statistics for ./po_97258_pan_0000000.tif 2021-01-19 17:54:05 {0} [ console ] : --> Reading statistics from file run_bundle/run-po_97258_pan_0000000-stats.tif 2021-01-19 17:54:05 {0} [ console ] : ./po_97258_pan_0000000.tif: [ lo: 0 hi: 255 mean: 73.9326 std_dev: 54.3488 ] 2021-01-19 17:54:05 {0} [ console ] : --> Computing statistics for ./po_97258_pan_0010000.tif 2021-01-19 17:54:05 {0} [ console ] : --> Reading statistics from file run_bundle/run-po_97258_pan_0010000-stats.tif 2021-01-19 17:54:05 {0} [ console ] : ./po_97258_pan_0010000.tif: [ lo: 0 hi: 255 mean: 72.0923 std_dev: 52.0436 ] 2021-01-19 17:54:05 {0} [ console ] : --> Matching interest points in StereoSession. 2021-01-19 17:54:05 {0} [ console ] : Using epipolar threshold = 188.656 2021-01-19 17:54:05 {0} [ console ] : IP uniqueness threshold = 0.8 2021-01-19 17:54:05 {0} [ console ] : Datum: Geodetic Datum --> Name: WGS_1984 Spheroid: WGS 84 Semi-major axis: 6378137 Semi-minor axis: 6356752.3142451793 Meridian: Greenwich at 0 Proj4 Str: +proj=longlat +datum=WGS84 +no_defs 2021-01-19 17:54:05 {0} [ console ] : Skipping rough homography. 2021-01-19 17:54:05 {0} [ console ] : Looking for IP in left image... 2021-01-19 17:54:05 {0} [ console ] : Using 1309 interest points per tile (1024^2 px). 2021-01-19 17:54:05 {0} [ console ] : Detecting IP 2021-01-19 17:54:06 {0} [ console ] : Building descriptors 2021-01-19 17:54:06 {0} [ console ] : Found interest points: 4999 2021-01-19 17:54:06 {0} [ console ] : Recording interest points to file: run_bundle/run-po_97258_pan_0000000.vwip 2021-01-19 17:54:06 {0} [ console ] : Looking for IP in right image... 2021-01-19 17:54:06 {0} [ console ] : Using 1309 interest points per tile (1024^2 px). 2021-01-19 17:54:06 {0} [ console ] : Detecting IP 2021-01-19 17:54:07 {0} [ console ] : Building descriptors 2021-01-19 17:54:07 {0} [ console ] : Found interest points: 4999 2021-01-19 17:54:07 {0} [ console ] : Recording interest points to file: run_bundle/run-po_97258_pan_0010000.vwip 2021-01-19 17:54:07 {0} [ console ] : --> Matching interest points using the epipolar line. 2021-01-19 17:54:07 {0} [ console ] : Uniqueness threshold: 0.8 2021-01-19 17:54:07 {0} [ console ] : Epipolar threshold: 188.656 2021-01-19 17:54:07 {0} [ console ] : Matching forward 2021-01-19 17:54:09 {0} [ console ] : ---> Obtained 4999 matches. 2021-01-19 17:54:09 {0} [ console ] : Matching backward 2021-01-19 17:54:10 {0} [ console ] : ---> Obtained 4999 matches. 2021-01-19 17:54:10 {0} [ console ] : Matched 1085 points. 2021-01-19 17:54:11 {0} [ console ] : Removed 31 points in stddev filtering. 2021-01-19 17:54:11 {0} [ console ] : Reduced matches to 1054 2021-01-19 17:54:11 {0} [ console ] : * Writing match file: run_bundle/run-po_97258_pan_0000000po_97258_pan_0010000.match 2021-01-19 17:54:11 {0} [ console ] : IP coverage fraction = 1 2021-01-19 17:54:11 {0} [ console ] : Number of matches in run_bundle/run-po_97258_pan_0000000__po_97258_pan_0010000.match 1054 2021-01-19 17:54:11 {0} [ console ] : Building the control network for run_bundle/run-po_97258_pan_0000000po_97258_pan_0010000.match. 2021-01-19 17:54:11 {0} [ console ] : Loading GCP files... 2021-01-19 17:54:11 {0} [ console ] : --> Bundle adjust pass: 0 2021-01-19 17:54:11 {0} [ console ] : Writing initial condition files... 2021-01-19 17:54:11 {0} [ console ] : Starting the Ceres optimizer... 2021-01-19 17:54:12 {0} [ console ] : Solver Summary (v 1.14.0-eigen-(3.3.7)-lapack-suitesparse-(5.6.0)-cxsparse-(3.2.0)-eigensparse-openmp-no_tbb)

                                 Original                  Reduced

Parameter blocks 1055 1055 Parameters 3171 3171 Residual blocks 2108 2108 Residuals 4224 4224

Minimizer TRUST_REGION

Dense linear algebra library EIGEN Trust region strategy LEVENBERG_MARQUARDT

                                    Given                     Used

Linear solver DENSE_SCHUR DENSE_SCHUR Threads 4 4 Linear solver ordering AUTOMATIC 1053,2 Schur structure 2,3,6 2,3,6

Cost: Initial 8.403343e+01 Final 2.009871e+01 Change 6.393472e+01

Minimizer iterations 17 Successful steps 12 Unsuccessful steps 5

Time (in seconds): Preprocessor 0.002329

Residual only evaluation 0.078715 (17) Jacobian & residual evaluation 1.097740 (12) Linear solver 0.068793 (17) Minimizer 1.254822

Postprocessor 0.000115 Total 1.257266

Termination: CONVERGENCE (Parameter tolerance reached. Relative step_norm: 1.887881e-09 <= 1.000000e-08.)

2021-01-19 17:54:12 {0} [ console ] : Writing final condition log files... 2021-01-19 17:54:12 {0} [ console ] : Removing pixel outliers in preparation for another solver attempt. 2021-01-19 17:54:12 {0} [ console ] : Outlier statistics: b = -0.133783, e = 0.249057. 2021-01-19 17:54:12 {0} [ console ] : Removing as outliers points with mean reprojection error > 1.5. 2021-01-19 17:54:12 {0} [ console ] : Removed 1 outliers by reprojection error. 2021-01-19 17:54:12 {0} [ console ] : Removed 0 outliers based on disparity of ip. 2021-01-19 17:54:12 {0} [ console ] : IP coverage fraction after cleaning = 1 2021-01-19 17:54:12 {0} [ console ] : Saving 1052 filtered interest points. 2021-01-19 17:54:12 {0} [ console ] : Writing: run_bundle/run-po_97258_pan_0000000__po_97258_pan_0010000-clean.match 2021-01-19 17:54:12 {0} [ console ] : --> Bundle adjust pass: 1 2021-01-19 17:54:12 {0} [ console ] : Starting the Ceres optimizer... 2021-01-19 17:54:14 {0} [ console ] : Solver Summary (v 1.14.0-eigen-(3.3.7)-lapack-suitesparse-(5.6.0)-cxsparse-(3.2.0)-eigensparse-openmp-no_tbb)

                                 Original                  Reduced

Parameter blocks 1054 1054 Parameters 3168 3168 Residual blocks 2106 2106 Residuals 4220 4220

Minimizer TRUST_REGION

Dense linear algebra library EIGEN Trust region strategy LEVENBERG_MARQUARDT

                                    Given                     Used

Linear solver DENSE_SCHUR DENSE_SCHUR Threads 4 4 Linear solver ordering AUTOMATIC 1052,2 Schur structure 2,3,6 2,3,6

Cost: Initial 8.324747e+01 Final 1.941424e+01 Change 6.383323e+01

Minimizer iterations 31 Successful steps 23 Unsuccessful steps 8

Time (in seconds): Preprocessor 0.001479

Residual only evaluation 0.154502 (31) Jacobian & residual evaluation 1.901398 (23) Linear solver 0.043280 (31) Minimizer 2.107783

Postprocessor 0.000081 Total 2.109343

Termination: CONVERGENCE (Parameter tolerance reached. Relative step_norm: 6.119454e-09 <= 1.000000e-08.)

2021-01-19 17:54:14 {0} [ console ] : Writing final condition log files... 2021-01-19 17:54:14 {0} [ console ] : Removing pixel outliers in preparation for another solver attempt. 2021-01-19 17:54:14 {0} [ console ] : Outlier statistics: b = -0.13281, e = 0.247542. 2021-01-19 17:54:14 {0} [ console ] : Removing as outliers points with mean reprojection error > 1.5. 2021-01-19 17:54:14 {0} [ console ] : Removed 0 outliers by reprojection error. 2021-01-19 17:54:14 {0} [ console ] : Removed 0 outliers based on disparity of ip. 2021-01-19 17:54:14 {0} [ console ] : IP coverage fraction after cleaning = 1 2021-01-19 17:54:14 {0} [ console ] : Saving 1052 filtered interest points. 2021-01-19 17:54:14 {0} [ console ] : Writing: run_bundle/run-po_97258_pan_0000000__po_97258_pan_0010000-clean.match 2021-01-19 17:54:14 {0} [ console ] : No new outliers removed, and the algorithm converged. No more passes are needed. 2021-01-19 17:54:14 {0} [ console ] : Writing: run_bundle/run-po_97258_pan_0000000.adjust 2021-01-19 17:54:14 {0} [ console ] : Writing: run_bundle/run-po_97258_pan_0010000.adjust

end

I have a more question about the outputs of the bundle adjustment ".adjust". I found there were only seven elements in the file, such as: "-1.88285418932678628 4.33993919949374529 -2.04130864424961533 0.999999993922814134 -6.07034995126527692e-06 4.56290377900619882e-06 -0.000109985009661008292 " There are 90 elements in the "_RPC.TXT" file. I am curious about that whether the above seven elements are the all parameters used to replace the 90 elements in the "*_RPC.txt" or there are other elements not showing.

adehecq commented 3 years ago

I did a side to side comparison of the log files (using meld) and the differences are very small. The number of IP matches is the same both for the initial set (1054) and after filtering (1052), so the problem does not seem to come from there. You can check that the match files are the same with stereo_gui: https://stereopipeline.readthedocs.io/en/latest/tools/stereo_gui.html?highlight=bundle_adjust#view-create-move-delete-save-interest-point-matches-gcp-and-vwip-files By the way, the third pass did not run as the algorithm converges at the second pass, so no need for num-passes 3. The solver also yield very similar output cost at both iterations. You can visually inspect the output of bundle_adjust by looking at the files initial/final_residuals_no_loss_function_pointmap_point_log.csv. You can convert them to raster files using ASP's point2dem with the following command: point2dem $outPrefix-initial_residuals_no_loss_function_pointmap_point_log.csv --csv-format '1:lon 2:lat 4:height_above_datum' --csv-proj4 '+proj=longlat +ellp s=WGS84 +datum=WGS84 +no_defs' (you can try setting your own output projection and resolution to improve the result, or you can open the files in a GIS)

Can you also show the difference of the .adjust files? And what is the range of the vertical offset you obtain between both results?

Oleg/Scott will have to answer about the content of the adjust file, I'm only using ASP with frame cameras...

zzqstar commented 3 years ago

Thank you for the reply. The values of the elements in two *.adjust files are: Adjustment results at first time: -0.740803432012757557 0.46522792383296413 -0.36206291138125779 0.999999997319504841 -6.6613907674501755e-06 -1.57762082126121799e-05 -7.1187972119639677e-05

Adjustment results at second time -1.88285418932678628 4.33993919949374529 -2.04130864424961533 0.999999993922814134 -6.07034995126527692e-06 4.56290377900619882e-06 -0.000109985009661008292

The offset of the two terrain models is about 2 meters at plane OXY. (The points were picked using Cloudcompare)

oleg-alexandrov commented 3 years ago

It is not surprising that two bundle adjustment passes would give different results. Bundle adjustment only makes the cameras correct relative to each other, not to the ground. Unless you use GCP. But even so, I thought bundle adjustment was deterministic. Maybe you can try running it with one thread only, so --num-threads 1 (if I recall right) and see if all runs give the same results.

oleg-alexandrov commented 3 years ago

Oh, and after bundle adjustment you better do some alignment to a known ground terrain, if you care about absolute location on the ground. We have the pc_align tool for that.

zzqstar commented 3 years ago

Thank you for the reply. I think that the comment I made about that the values of the elements in two *.adjust files were different

Thank you for the reply. The values of the elements in two *.adjust files are: Adjustment results at first time: -0.740803432012757557 0.46522792383296413 -0.36206291138125779 0.999999997319504841 -6.6613907674501755e-06 -1.57762082126121799e-05 -7.1187972119639677e-05

Adjustment results at second time -1.88285418932678628 4.33993919949374529 -2.04130864424961533 0.999999993922814134 -6.07034995126527692e-06 4.56290377900619882e-06 -0.000109985009661008292

The offset of the two terrain models is about 2 meters at plane OXY. (The points were picked using Cloudcompare)

may cause misunderstanding. To be clear, the above adjustment results at first time and second time are not related to the paramter --num-passes in bundle adjustment process, and not the results corresponding to the --num-passes set to be 1 and 2, seperately. The -num-passes were set to be 1, and bundle adjustment were excuted two times seperately (completely indepentent), generating two different results. To be more clealy, the bundle adjustment were excuted seven times (completely independent), followed by stereo matching process, generating seven results. The centroids of the seven point clouds were estimated as follows (the centroid may not be a very rigorous index for describing the offset among these results): X Y Z 485953.308387414 3620253.28706555 12.3252335959041 485952.967302522 3620252.88273437 10.1524395843361 485953.169409302 3620249.52577605 4.20331607201715 485956.01107743 3620247.27427926 14.3982800668468 485954.023897953 3620250.62003086 9.25236287007571 485956.153328523 3620251.93539036 12.8895646558363 485954.741403963 3620252.45207816 12.3219433726963

I have a more question about the outputs of the bundle adjustment ".adjust". I found there were only seven elements in the file, such as: "-1.88285418932678628 4.33993919949374529 -2.04130864424961533 0.999999993922814134 -6.07034995126527692e-06 4.56290377900619882e-06 -0.000109985009661008292 " There are 90 elements in the "_RPC.TXT" file. I am curious about that whether the above seven elements are the all parameters used to replace the 90 elements in the "*_RPC.txt" or there are other elements not showing.

zzqstar commented 3 years ago

Thank you very much, Oleg. The parameter --threads should be set to be 1 in bundle adjustment, then, the bundle adjustment is deterministic. Default value for --threads will caused an uncertainty of the result (I guess that multithread process is not suitable the removal of the false matching point pairs, if RANSAC was used ) .

oleg-alexandrov commented 3 years ago

Glad to hear that --threads 1 makes bundle adjustment deterministic. We kind of knew that, but I forgot before.

The .adjust file has a translation vector and a rotation encoded as a quaternion. We do not float individually the 90 elements in the RPC file. In fact we keep them all fixed, and rigidly rotate and translate the camera determined by these parameters.

adehecq commented 3 years ago

What is the reason behind the use of multi-threading causing an undeterministic behavior?

ScottMcMichael commented 3 years ago

This is a limitation of the Ceres solver we use.

Scott


From: Amaury Dehecq notifications@github.com Sent: Tuesday, January 19, 2021 11:57 PM To: NeoGeographyToolkit/StereoPipeline StereoPipeline@noreply.github.com Cc: ScottMcMichael SMcMichael@sgt-inc.com; Comment comment@noreply.github.com Subject: [External] Re: [NeoGeographyToolkit/StereoPipeline] loading RPC model when stereo with Map-projected Images v2.6.1 (#221)

What is the reason behind the use of multi-threading causing an undeterministic behavior?

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fprotect2.fireeye.com%2Fv1%2Furl%3Fk%3D3eacb9b5-613780c2-3eacc819-8620f0aedd21-976682b1f507197a%26q%3D1%26e%3Dca9a9c2c-5d24-46b2-8216-aa95ccfaf72d%26u%3Dhttps%253A%252F%252Fgithub.com%252FNeoGeographyToolkit%252FStereoPipeline%252Fissues%252F221%2523issuecomment-763411642&data=04%7C01%7Cscott.t.mcmichael%40nasa.gov%7C73bcdd0e94cc488b335108d8bd191d2a%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C637467262860336545%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=Rv9fIqGe7PFn87FaF6TzxE9Rv%2FnzdMi002lVk0GnMnU%3D&reserved=0, or unsubscribehttps://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fprotect2.fireeye.com%2Fv1%2Furl%3Fk%3D65e28e29-3a79b75e-65e2ff85-8620f0aedd21-cac539c2c5cfe2cd%26q%3D1%26e%3Dca9a9c2c-5d24-46b2-8216-aa95ccfaf72d%26u%3Dhttps%253A%252F%252Fgithub.com%252Fnotifications%252Funsubscribe-auth%252FABDVLRD7SAAE2NNRKHXRX5DS22EFPANCNFSM4GRHMLEA&data=04%7C01%7Cscott.t.mcmichael%40nasa.gov%7C73bcdd0e94cc488b335108d8bd191d2a%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C637467262860346510%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=GViz2d%2FaGEZe2YPRQj8qswwEG%2B5CptOmi91XF0wLHAQ%3D&reserved=0.

This e-mail, including any attached files, may contain confidential information, privileged information and/or trade secrets for the sole use of the intended recipient. Any review, use, distribution, or disclosure by others is strictly prohibited. If you are not the intended recipient (or authorized to receive information for the intended recipient), please contact the sender by reply e-mail and delete all copies of this message.

zzqstar commented 3 years ago

Thank you for the reply. There is another question (report a error) while using the parallel_stereo commands in ASP to stereo matching a satellite image (the number of pixels is about 35000*40000). The stereo.default is as follows, stereo.default: stereo-algorithm 1 cost-mode 4 corr-kernel 7 7 subpixel-mode 12 median-filter-size 3 alignment-method affineepipolar

Error Information: Error This error is located in /AmesStereoPipeline/libexec/stereoutils.py, at line 114: # TODO: Move this to asp_system_utils

A very simple wrapper around subprocess

def generic_run(cmd, verbose): cmd_str = asp_string_utils.argListToString(cmd) if verbose: print(cmd_str)

try: code = subprocess.call(cmd) #line 110 except OSError as e: raise Exception('%s: %s' % (cmd_str, e)) if code != 0: #code != 0: raise Exception('Failed to run: ' + cmdstr) #Line 114 The code is not not equal to zero. This error appeared after the Stage 1 -- Correlation Finished, and before the Stage 2 -- Blending. I have not capable of debugging the program at line 110 (Not familar with python, do most job using C++ and matlab...) currently. The problem is about paramter "cmd" at line 110. However, the above error will not appear, if using the stereo.default as follows stereo.default: stereo-algorithm 0 cost-mode 2 corr-kernel 21 21 subpixel-mode 2 median-filter-size 3 alignment-method affineepipolar

oleg-alexandrov commented 3 years ago

Something failed during stereo correlation and your screenshot is not enough to tell what it was.

Did it happen with any other dataset? You can also debug this on a small clip, by using gdal_translate to make clips from both the left and the right images while keeping the (0, 0) upper-left corner in both, and trying parallel_stereo with these clips.

This code did not change for a very long time, so it is likely a user error. Maybe one of your nodes ran out of memory and died. You will need to study this carefully to understand the point of failure.

ScottMcMichael commented 3 years ago

SGM can use a lot of memory so with 16 processes there is a good chance you ran out of memory. You can increase the log level in your ~/.vwrc file to see more messages about how memory is being consumed, and you should be able to find some threads discussing the issue in our support forum: https://groups.google.com/g/ames-stereo-pipeline-support

Scott


From: Oleg Alexandrov notifications@github.com Sent: Friday, January 22, 2021 9:39 AM To: NeoGeographyToolkit/StereoPipeline StereoPipeline@noreply.github.com Cc: ScottMcMichael SMcMichael@sgt-inc.com; Comment comment@noreply.github.com Subject: [External] Re: [NeoGeographyToolkit/StereoPipeline] loading RPC model when stereo with Map-projected Images v2.6.1 (#221)

Something failed during stereo correlation and your screenshot is not enough to tell what it was.

Did it happen with any other dataset? You can also debug this on a small clip, by using gdal_translate to make clips from both the left and the right images while keeping the (0, 0) upper-left corner in both, and trying parallel_stereo with these clips.

This code did not change for a very long time, so it is likely a user error. Maybe one of your nodes ran out of memory and died. You will need to study this carefully to understand the point of failure.

— You are receiving this because you commented. Reply to this email directly, view it on GitHubhttps://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fprotect2.fireeye.com%2Fv1%2Furl%3Fk%3D5fe2d3de-0079eb14-5fe2a272-86324fd7403a-f0d8605946100d13%26q%3D1%26e%3D9f238024-1f02-4496-9a5e-ed03c6a670fd%26u%3Dhttps%253A%252F%252Fgithub.com%252FNeoGeographyToolkit%252FStereoPipeline%252Fissues%252F221%2523issuecomment-765575805&data=04%7C01%7Cscott.t.mcmichael%40nasa.gov%7C735c9f6b999a426fe3f108d8befccbf1%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C637469340254055186%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=um%2BkIvdV%2Fej0vqc6t0HfplSrlEIv0Jfmer6WmlZEaxI%3D&reserved=0, or unsubscribehttps://gcc02.safelinks.protection.outlook.com/?url=https%3A%2F%2Fprotect2.fireeye.com%2Fv1%2Furl%3Fk%3D22fff25b-7d64ca91-22ff83f7-86324fd7403a-ffbd809100675c2d%26q%3D1%26e%3D9f238024-1f02-4496-9a5e-ed03c6a670fd%26u%3Dhttps%253A%252F%252Fgithub.com%252Fnotifications%252Funsubscribe-auth%252FABDVLRHNH2A3TP4DNJMMSU3S3GZ3RANCNFSM4GRHMLEA&data=04%7C01%7Cscott.t.mcmichael%40nasa.gov%7C735c9f6b999a426fe3f108d8befccbf1%7C7005d45845be48ae8140d43da96dd17b%7C0%7C0%7C637469340254055186%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=CZ2dwAGsdZ2ABGQ%2Bykyj%2FtSmGIxl9zwiMHePHtiw64c%3D&reserved=0.

This e-mail, including any attached files, may contain confidential information, privileged information and/or trade secrets for the sole use of the intended recipient. Any review, use, distribution, or disclosure by others is strictly prohibited. If you are not the intended recipient (or authorized to receive information for the intended recipient), please contact the sender by reply e-mail and delete all copies of this message.

zzqstar commented 3 years ago

Thank you. The problem with stereo correlation was solved by lowering the value of prameter --processes. There are 32G memory volume of my computer. To make sure that memory exceptions will never happen, the --processes is set to be 4. The memory consumtion was also monitored, and none exception happened. However, there is an error report at the very beginning of the stage 3 ---Filtering (After the blending): GdalIO: '*.tif' not recognised as a supported file format (code =4) Error-Filtering I started a new norminal window, and used the "gdalinfo" command to read the infomation of the image not recognized above, no error reported. Using a smaller and different image, it did not report a problem likewise.

zzqstar commented 3 years ago

Glad to hear that --threads 1 makes bundle adjustment deterministic. We kind of knew that, but I forgot before.

The .adjust file has a translation vector and a rotation encoded as a quaternion. We do not float individually the 90 elements in the RPC file. In fact we keep them all fixed, and rigidly rotate and translate the camera determined by these parameters.

Dear Olgy, sorry to disturb you again. The .adjust files consisting of seven elements were generated using the function bundle_adjustment, and seven elements (transformation vector and a quaternion) can be used to compensate the raw rational polynomial coefficients (RPC). I am a little confused about how to use these seven elements to calculate the pixel coordinates cooresponding to the ground position. QQ截图20210517121449 Or, if the .adjust + old_RPC.txt can be converted into the new RPC file (new_RPC.txt).

oleg-alexandrov commented 3 years ago

We don't use the bundle adjustment values to modify these ratios. The logic we use is the following:

Given a point xyz in space that we want to project into the camera, we first apply a small rotation and translation to that xyz using the 7 coefficients in the .adjust file (4 for rotation and 3 for translation). The we plug the new xyz into the above formulas for ln and sn.

zzqstar commented 3 years ago

Dear Olgy, thank you very much for the reply. As you have said, I computed the pixel coordinates corresponding to the ground position using the bundle adjustment results and the raw RPC parameters. However, there are some problem:

Given a ground point (X0, Y0, Z0) (ground projected coordinates), and its corrsponding pixel coordinates true value in image are (Line0, Sample0) (no mor thant 1 pixel error). Using the RPC parameters only to compute the pixel coordinates (Line1, Sample1) of ground point (X0, Y0, Z0), the difference between (Line0, Sample0) and (Line1, Sample1) will be within 10 pixels.

While using the RPC parameters and bundle adjustment results (a small rotation and translation applied to XYZ before the above ratio formulas), the difference between the final pixel coordinates (Line2, Sample2) and (Line0, Sample0) will be hundreds to thousands pixels. But, based on the the bundle adjustment results (compared with none bundle adjustment), the higher accuracy of the DEM (generated with the commands parallel_stereo and point2dem) has been achieved. Thus, I doubt if there were some mistakes to compute the (Line2, Sample2):

"-1.88285418932678628 4.33993919949374529 -2.04130864424961533 0.999999993922814134 -6.07034995126527692e-06 4.56290377900619882e-06 -0.000109985009661008292 " The three elements in the first line are the translation vector (V). The second line is a quaternion (q0, q1, q2, q3), and can be coverted into a roataion matrix (R). The transfomation applied to ground point can be: R[X0 Y0 Z0]'+V, R[X0 Y0 Z0]'-V, R([X0 Y0 Z0]'+V), or R([X0 Y0 Z0]'-V), and result can be coverted into latitude, longitude and height before using above rpc ratio formulas. I tried these four transformation formulas, and the final results of pixel coordinates were all deviating far from the true value (Line0, Sample0).

oleg-alexandrov commented 3 years ago

There is probably some issue with the input data or how you do things, I believe. All this logic has been there for a long time and is rather well tested.

A simple solution out of the whole thing is to not use GCP. I would just use bundle adjustment (say with a larger --camera-weight if you notice that bundle adjustment moves the cameras too much), then create a DEM, and compare that one to a different DEM of the area, such as SRTM. Then one can use pc_align to bring the ASP DEM to the desired location.

If you want to use GCP, I wonder if you used our stereo_gui tool for creating GCP.

You can also start practicing by creating GCP not based on prior knowledge, but using stereo_gui with an orthoimage obtained by mapprojecting an image onto a DEM ASP created to start with. This way you'll know the answer, the GCP should make little difference as the ASP DEM is already self-consistent with the cameras, if not consistent with the true position on the ground.

Lastly, I wonder if you looked at the no_loss_pointmap output file created by bundle adjustment,. It shows the residuals.

Hopefully after you practice more the results may improve. Or else I would need a detailed log of exactly what you tried to do and the results you got to figure out what is going on.

zzqstar commented 3 years ago

Thank you very much for the reply.

I tried to use the bundle adjustment results (.adjust). Briefly, given a ground point (lat,lon, h), I convert this point into cartesian coordinate (xyz), and then apply this point (xyz) with a rotation and translation (seven elements in .adjust), and finally, plug the new xyz (will be converted into ellipsoidal coordinates) into the rpc ratio formulas to compute the pixel coordinate (Is this the right method? ). The true value for the corrsponding pixel coordinate is deterministic, and the calculation result deviate far from the true value (hundreds to thousands pixels). The difference will be less than 20 pixels between the true value with the calculation result if I plug the (lat, lon, h) directly into rpc ratio formulas (don't use the bundle adjustment results).

Following are the bundle adjustment logs: _/home/zzqstar/Desktop/StereoPipeline270/libexec/bundle_adjust ./BWD.tif ./FWD.tif ground_control_points.gcp --datum WGS_1984 --num-passes 2 --ip-num-ransac-iterations 1500 -r 30 --remove-outliers-params 75 2.5 1.5 2.5 -o run_bundle/run --fix-gcp-xyz Vision Workbench log started at 2021-05-20 19:12:30. 2021-05-20 19:12:30 {0} [ console ] : Using session: rpc. 2021-05-20 19:12:30 {0} [ console ] : --> Computing statistics for ./BWD.tif 2021-05-20 19:12:30 {0} [ console ] : using downsample scale: 38 2021-05-20 19:12:31 {0} [ console ] : Writing stats file: run_bundle/run-BWD-stats.tif 2021-05-20 19:12:31 {0} [ console ] : ./BWD.tif 2021-05-20 19:12:31 {0} [ console ] : Using session: rpc. 2021-05-20 19:12:31 {0} [ console ] : --> Computing statistics for ./FWD.tif 2021-05-20 19:12:31 {0} [ console ] : using downsample scale: 32 2021-05-20 19:12:31 {0} [ console ] : Writing stats file: run_bundle/run-FWD-stats.tif 2021-05-20 19:12:31 {0} [ console ] : ./FWD.tif 2021-05-20 19:12:31 {0} [ console ] : Using session: rpc. 2021-05-20 19:12:31 {0} [ console ] : Loading camera model: ./BWD.tif 2021-05-20 19:12:31 {0} [ console ] : Using session: rpc. 2021-05-20 19:12:31 {0} [ console ] : Loading camera model: ./FWD.tif 2021-05-20 19:12:31 {0} [ console ] : Using session: rpc. 2021-05-20 19:12:31 {0} [ console ] : Using session: rpc. 2021-05-20 19:12:31 {0} [ console ] : --> Computing statistics for ./BWD.tif 2021-05-20 19:12:31 {0} [ console ] : --> Reading statistics from file run_bundle/run-BWD-stats.tif 2021-05-20 19:12:31 {0} [ console ] : ./BWD.tif 2021-05-20 19:12:31 {0} [ console ] : --> Computing statistics for ./FWD.tif 2021-05-20 19:12:31 {0} [ console ] : --> Reading statistics from file run_bundle/run-FWD-stats.tif 2021-05-20 19:12:31 {0} [ console ] : ./FWD.tif 2021-05-20 19:12:31 {0} [ console ] : --> Matching interest points in StereoSession. 2021-05-20 19:12:31 {0} [ console ] : Warning: The RPC model minimum height is above the zero datum. 2021-05-20 19:12:31 {0} [ console ] : RPC model min and max heights above datum: 1243.18 1646.4. 2021-05-20 19:12:31 {0} [ console ] : Adjusting the datum to compensate, for the purpose of alignment. 2021-05-20 19:12:31 {0} [ console ] : The new datum height will be at 1444.79 relative to the previous one. 2021-05-20 19:12:31 {0} [ console ] : Using epipolar threshold = 3581.57 2021-05-20 19:12:31 {0} [ console ] : IP uniqueness threshold = 0.8 2021-05-20 19:12:31 {0} [ console ] : Datum: Geodetic Datum --> Name: WGS_1984 Spheroid: WGS 84 Semi-major axis: 6379581.79 Semi-minor axis: 6358197.1042451793 Meridian: Greenwich at 0 Proj4 Str: +a=6.37958e+06 +b=6.3582e+06 2021-05-20 19:12:31 {0} [ console ] : Skipping rough homography. 2021-05-20 19:12:31 {0} [ console ] : Looking for IP in left image... 2021-05-20 19:12:31 {0} [ console ] : Using 50 interest points per tile (1024^2 px). 2021-05-20 19:12:31 {0} [ console ] : Detecting IP 2021-05-20 19:16:42 {0} [ console ] : Building descriptors 2021-05-20 19:16:54 {0} [ console ] : Found interest points: 68469 2021-05-20 19:16:54 {0} [ console ] : Recording interest points to file: run_bundle/run-BWD.vwip 2021-05-20 19:16:54 {0} [ console ] : Looking for IP in right image... 2021-05-20 19:16:54 {0} [ console ] : Using 50 interest points per tile (1024^2 px). 2021-05-20 19:16:54 {0} [ console ] : Detecting IP 2021-05-20 19:19:56 {0} [ console ] : Building descriptors 2021-05-20 19:20:05 {0} [ console ] : Found interest points: 46238 2021-05-20 19:20:05 {0} [ console ] : Recording interest points to file: run_bundle/run-FWD.vwip 2021-05-20 19:20:05 {0} [ console ] : --> Matching interest points using the epipolar line. 2021-05-20 19:20:05 {0} [ console ] : Uniqueness threshold: 0.8 2021-05-20 19:20:05 {0} [ console ] : Epipolar threshold: 3581.57 2021-05-20 19:20:05 {0} [ console ] : Matching forward 2021-05-20 19:22:39 {0} [ console ] : ---> Obtained 68469 matches. 2021-05-20 19:22:39 {0} [ console ] : Matching backward 2021-05-20 19:25:36 {0} [ console ] : ---> Obtained 46238 matches. 2021-05-20 19:25:36 {0} [ console ] : Matched 6083 points. 2021-05-20 19:25:53 {0} [ console ] : Removed 231 points in stddev filtering. 2021-05-20 19:25:53 {0} [ console ] : Reduced matches to 5852 2021-05-20 19:25:53 {0} [ console ] : * Writing match file: run_bundle/run-BWDFWD.match 2021-05-20 19:25:53 {0} [ console ] : IP coverage fraction = 0.568498 2021-05-20 19:25:53 {0} [ console ] : Number of matches in run_bundle/run-BWD__FWD.match 5852 2021-05-20 19:25:53 {0} [ console ] : Building the control network for run_bundle/run-BWDFWD.match. 2021-05-20 19:25:54 {0} [ ba ] : Warning: 6 control points removed due to spiral errors. 2021-05-20 19:25:54 {0} [ console ] : Loading GCP files... 2021-05-20 19:25:54 {0} [ console ] : Loading GCP file: ground_control_points.gcp 2021-05-20 19:25:54 {0} [ console ] : WARNING: GCPs are over 100 km from the other points. Are your lat/lon GCP coordinates swapped? 2021-05-20 19:25:54 {0} [ console ] : --> Bundle adjust pass: 0 2021-05-20 19:25:54 {0} [ console ] : Writing initial condition files... 2021-05-20 19:25:54 {0} [ console ] : Starting the Ceres optimizer... 2021-05-20 19:26:03 {0} [ console ] : Solver Summary (v 1.14.0-eigen-(3.3.7)-lapack-suitesparse-(5.6.0)-cxsparse-(3.2.0)-eigensparse-openmp-no_tbb)

                             Original                  Reduced

Parameter blocks 5858 5842 Parameters 17580 17532 Residual blocks 11730 11714 Residuals 23484 23436

Minimizer TRUST_REGION

Dense linear algebra library EIGEN Trust region strategy LEVENBERG_MARQUARDT

                                Given                     Used

Linear solver DENSE_SCHUR DENSE_SCHUR Threads 1 1 Linear solver ordering AUTOMATIC 5840,2 Schur structure 2,3,6 2,3,6

Cost: Initial 4.546920e+03 Final 2.928999e+02 Change 4.254020e+03

Minimizer iterations 22 Successful steps 22 Unsuccessful steps 0

Time (in seconds): Preprocessor 0.009602

Residual only evaluation 0.492364 (22) Jacobian & residual evaluation 8.733223 (22) Linear solver 0.111321 (22) Minimizer 9.365918

Postprocessor 0.000540 Total 9.376060

Termination: CONVERGENCE (Parameter tolerance reached. Relative step_norm: 5.569133e-09 <= 1.000000e-08.)

2021-05-20 19:26:03 {0} [ console ] : Writing final condition log files... 2021-05-20 19:26:03 {0} [ console ] : Ground control point results: 2021-05-20 19:26:03 {0} [ console ] : input_gcp optimized_gcp diff 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : xyz: 2021-05-20 19:26:03 {0} [ console ] : llh: 2021-05-20 19:26:03 {0} [ console ] : Removing pixel outliers in preparation for another solver attempt. 2021-05-20 19:26:03 {0} [ console ] : Outlier statistics: b = -0.151751, e = 0.268258. 2021-05-20 19:26:03 {0} [ console ] : Removing as outliers points with mean reprojection error > 1.5. 2021-05-20 19:26:03 {0} [ console ] : Removed 97 outliers by reprojection error. 2021-05-20 19:26:03 {0} [ console ] : Removed 0 outliers based on disparity of ip. 2021-05-20 19:26:03 {0} [ console ] : IP coverage fraction after cleaning = 0.834444 2021-05-20 19:26:03 {0} [ console ] : Saving 5743 filtered interest points. 2021-05-20 19:26:03 {0} [ console ] : Writing: run_bundle/run-BWD__FWD-clean.match 2021-05-20 19:26:03 {0} [ console ] : --> Bundle adjust pass: 1 2021-05-20 19:26:04 {0} [ console ] : Starting the Ceres optimizer... 2021-05-20 19:26:08 {0} [ console ] : Solver Summary (v 1.14.0-eigen-(3.3.7)-lapack-suitesparse-(5.6.0)-cxsparse-(3.2.0)-eigensparse-openmp-no_tbb)

                             Original                  Reduced

Parameter blocks 5761 5745 Parameters 17289 17241 Residual blocks 11536 11520 Residuals 23096 23048

Minimizer TRUST_REGION

Dense linear algebra library EIGEN Trust region strategy LEVENBERG_MARQUARDT

                                Given                     Used

Linear solver DENSE_SCHUR DENSE_SCHUR Threads 1 1 Linear solver ordering AUTOMATIC 5743,2 Schur structure 2,3,6 2,3,6

Cost: Initial 4.311028e+03 Final 1.599738e+02 Change 4.151054e+03

Minimizer iterations 11 Successful steps 11 Unsuccessful steps 0

Time (in seconds): Preprocessor 0.008426

Residual only evaluation 0.240803 (11) Jacobian & residual evaluation 4.265211 (11) Linear solver 0.054686 (11) Minimizer 4.574986

Postprocessor 0.000513 Total 4.583925

Termination: CONVERGENCE (Parameter tolerance reached. Relative step_norm: 9.236386e-09 <= 1.000000e-08.)

2021-05-20 19:26:08 {0} [ console ] : Writing final condition log files... 2021-05-20 19:26:08 {0} [ console ] : Ground control point results: 2021-05-20 19:26:08 {0} [ console ] : input_gcp optimized_gcp diff 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : xyz: 2021-05-20 19:26:08 {0} [ console ] : llh: 2021-05-20 19:26:08 {0} [ console ] : Removed 0 outliers based on disparity of ip. 2021-05-20 19:26:08 {0} [ console ] : IP coverage fraction after cleaning = 0.834444 2021-05-20 19:26:08 {0} [ console ] : Saving 5743 filtered interest points. 2021-05-20 19:26:08 {0} [ console ] : Writing: run_bundle/run-BWD__FWD-clean.match 2021-05-20 19:26:08 {0} [ console ] : Writing: run_bundle/run-BWD.adjust 2021-05-20 19:26:08 {0} [ console ] : Writing: runbundle/run-FWD.adjust

The GCP file was add in the above bundle adjustment (Even withour GCP involved in bundle adjustment, the above calculation results (pixel coordinates) still deviate far from the true value). The stereo_gui tool was not used for creating GCP. however, the points in GCP file were strictly arranged according to the asp book. The residuals for the points (including the added gcp point ) in final_residuals_no _less_function_pointmap are relatively small (less than 5).

oleg-alexandrov commented 3 years ago

I am not sure I understand what you are doing and what you mean by "results (pixel coordinates) still deviate far from the true value".

As before, I will suggest avoiding GCP altogether, run bundle adjustment, run stereo, and compare your obtained terrain to any previous terrain you have. Also run stereo with and without bundle adjustment, and look at the point2dem output file having the intersection error image (see point2dem --errorimage). Bundle adjustment is meant to make the cameras more consistent.

Briefly, given a ground point (lat,lon, h), I convert this point into cartesian coordinate (xyz), and then apply this point (xyz) with a rotation and translation (seven elements in *.adjust), and finally, plug the new xyz (will be converted into ellipsoidal coordinates) into the rpc ratio formulas to compute the pixel coordinate (Is this the right method? ).

It is a bit more hairy than that, because you don't know around which point the rotation happens. That is something obscure in our code which you don't have access to, it is the camera center at pixel (0, 0), which for RPC models is computed in some arbitrary way.

So what I said about .adjust files should be interpreted in the most general way.

zzqstar commented 3 years ago

Dear Oleg, thank you for the reply. I may have got what you have said. I want to calculate the image pixel coordinates by using the ground coordinates (xyz), RPC models, and *.adjustment exports. However, this may be infeasible for me because I don't have the access to the rotation center for ground coordinates (xyz), as you have said. I take the rotation center as (0,0,0) in the previous calculation, which lead a large discrepancy between the calculated pixel coordinates and its exact (true) location.

'That is something obscure in our code which you don't have access to, it is the camera center at pixel (0, 0), which for RPC models is computed in some arbitrary way.' I not very sure if the arbitary way for RPC models to compute camera center at pixel (0,0) in bundle adjustment will cause some errors. In some cases, the residuals for the GCP in final_residuals_no _less_function_pointmap are not very small (up to 30), and the precision for the final DSM is even decreased, compared with that of stereo matching generation without GCP.

oleg-alexandrov commented 3 years ago

If you are getting big residuals in the final_residuals_no _less_function_pointmap output file (like 30), that likely means the GCP are not good, or inconsistent with each other. You can always try to increase the --robust-threshold, say to 5, to force bundle adjustment to optimize harder these residuals at the expense of other things, but that is more of an artificial solution.

You can also set --camera-weight 0 so that the cameras are allowed to move as much as they want in order to optimize your GCP. That is also a bit of a dangerous thing to do.

As before, your best bet is to first run bundle adjustment with no gcp, run stereo, run point2dem --errorimage, then do gdalinfo -stats on the obtained error image file. That error is the minimal distance between rays. That error image can be also colorized, such as colormap --min -0 --max 5 error.tif. This is useful in understanding if your RPC images and cameras are at least self-consistent. That is a very important sanity check.

Then, if you get a DEM which looks good but is misaligned, doing pc_align is helpful, followed by using geodiff --absolute to find the absolute difference between the two DEMs. The pc_align tool is better than GCP because it uses a lot of dense points and can remove outliers, while GCP are very sparse.

If all these give great results, but GCP does not, that will mean your GCP are not so good.

zzqstar commented 3 years ago

Thank you vey much, Oldg. After many repeated tests, I have finally found that the problem occurs in the GCP for whcich the residuals in final_residuals_no _less_function_pointmap are not very small (like 30) in previous tests. The bundle_adjust works well after the update of the GCP.

You mentioned that the camera center at pixel (0, 0) in bundle adjustment process, if for RPC models, is computed in some arbitrary way. Can I ask you how this arbitrary way was? And how to calculate the camer center as pixel (0,0 ) for RPC models? I tried to calculate this rotation center by using pixel(0,0) + RPC parameters + height_offset (in PRC models, assume to be the height at pixel (0,0)). However, taking the height_offset parameter as the height for pixel(0,0) seems to be improper. The image pixel coordinates calculated by using the ground coordinates (xyz), RPC models, and *.adjustment exports (the rotation center was calculated using the above method) are still deviating away from its true position (about 20 pixels). But this deviation has been much smaller than that of the previous calculated results (hundreds to thousands pixels, without calculation of rotation center).

Thank you again for the previous reply and the help.

oleg-alexandrov commented 3 years ago

Here's the code for the formula for the camera center for RPC cameras.

https://github.com/NeoGeographyToolkit/StereoPipeline/blob/master/src/asp/Camera/RPCModel.cc#L553

You'll need to trace the logic through this piece of code, which will require some C++ knowledge. The calculation is rather straightforward though, it uses the formulas for RPC coefficients, and the fact that most of the logic happens in the datum-based coordinate system, so where instead of x, y, z relative to ECEF (earth center) the values of lon, lan, and height above datum are used.

ShashankBice commented 3 years ago

Hi folks, One trick which might work to get a sane RPC model with the adjust term incorporated is to use the cam2rpc program, which can accept the --bundle-adjust-prefix. This should produce a new set of RPC model which should be good enough.. The caveat here is that the program crops the images a bit, so you should take note of the image dimensions which were used to compute the RPC model during cam2rpc invocation and then crop your original image to that extent using the gdal_translate command. At this point, you should be getting the world to image transformations correctly for your GCPs using the cam2rpc derived RPC, without worrying about the adjust files. That will be when you get the bundle adjustment right though ;) If my long sentences were a bit hard to follow, I am sure Oleg can explain this in a better way :)

Shashank

On Tue, May 25, 2021 at 9:16 AM Oleg Alexandrov @.***> wrote:

Here's the code for the formula for the camera center for RPC cameras.

https://github.com/NeoGeographyToolkit/StereoPipeline/blob/master/src/asp/Camera/RPCModel.cc#L553

You'll need to trace the logic through this piece of code, which will require some C++ knowledge. The calculation is rather straightforward though, it uses the formulas for RPC coefficients, and the fact that most of the logic happens in the datum-based coordinate system, so where instead of x, y, z relative to ECEF (earth center) the values of lon, lan, and height above datum are used.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/NeoGeographyToolkit/StereoPipeline/issues/221#issuecomment-848013231, or unsubscribe https://github.com/notifications/unsubscribe-auth/AG5K5UV3TSRO7Q3Q5HT45S3TPPEFPANCNFSM4GRHMLEA .

-- Shashank Bhushan Graduate Student/Research Assistant Civil & Environmental Engineering University of Washington, Seattle

ShashankBice commented 3 years ago

We can skip the gdal_translate step for cropping images, by using this flag --save-tif-image in the cam2rpc command itself.

zzqstar commented 3 years ago

Thank you very much, Oleg. The problem has been solved after an analysis of the code as you suggested.

zzqstar commented 3 years ago

Thank you for your reply, ShahankBice. I also tried to use the tool cam2rpc. However, this commands needs the range (longitude, latitude and height) of the surface covered by the satellite image. The generated RPC parameters will be less accuracte if the range values are not precise enough.

oleg-alexandrov commented 3 years ago

Actually the RPC approximation should do rather well if your estimated lon-lat-height range box is larger than the true box, say by a factor of 2, because RPC uses high-degree polynomials so it should be able to fit well a bigger area. Though of course if it becomes too big the accuracy will suffer.

Glad the code was useful.