Pycuampcor v2 #922

lijun99 · 2025-01-16T23:40:30Z

PyCuAmpcor Version 2 introduces a new workflow inspired by GrIMP/SpeckleSource, in addition to the original ROIPAC/ampcor workflow. Below is an overview of the new workflow and additional improvements:

Key Features of the New GrIMP/SpeckleSource Workflow:

Starts with an anti-aliasing oversampled secondary chip over the entire search range.
Adds extra margins to facilitate correlation surface extraction.
Offers optional double-precision computations (also extended to the ROIPAC/ampcor workflow for consistency).

The GrIMP/SpeckleSource workflow may be slower than ROIPAC/ampcor due to increased computational demands.
However, it may improve accuracy in scenarios with low signal-to-noise ratio (SNR). More details on the workflows and how to choose can be found at contrib/PyCuAmpcor/README.md.

Other Improvements

Introduced an option to fix the starting pixel location, which is ideal for comparing results from different window sizes.
Added an example script for plotting the offset field, making it easier for users to visualize results.
Included correlation peak values in the output, alongside SNR and covariance, providing additional accuracy measures.
Resolved a memory deallocation issue that occurred when the run procedure was called repeatedly within a single script.

in this update, * either the starting pixel(corner) or the center of window location may be specified (see README) * by specifying the starting/ending pixel locations, the reference/search window may go beyond the range of the original images. Zeros are filled in instead. In this case, the offset fields at the edge may be incorrect. * The range check is no longer needed. It only runs in DEBUG mode and issues warning messages when pixels are outside the image range.

…flag

…han 640 with DP

…memory

…tional approach to side Nyquist frequency (N/2) with negative frequency

…back

…virtual class for different workflows

…ing back Sinc Oversampler

hfattahi · 2025-01-17T05:42:00Z

Thank you @lijun99 for this extensive PR. Based on our offline discussions I'm aware of most of developments in this PR. However, for clarity and to help the user community following this repository would you please update the description of the PR and provide some explanation in few bullets summarizing the changes that comes in this PR.

vbrancat

@lijun99 Thank you for this PR! It’s very exciting, and I’m eager to test it. I will leave some further comments after I am done with the testing.

In the meantime, I’ve added some comments and questions on the code. Since I’m not a CUDA programmer, I’ll leave it to @gmgunter and @rtburns-jpl to weigh in on the CUDA-specific parts. Overall, the PR looks great to me, but there are a few comments and considerations that I think need to be addressed:

I’m not a fan of the ROIPAC vs. GRiMP terminology, as I believe most users won’t be familiar with it either. At a high level (e.g., in the Python interface), we should frame this in terms of enabling or disabling anti-aliasing oversampling and provide a clear explanation of the advantages and disadvantages of each option. This approach would ensure complete transparency for users.
It’s unclear to me how offset estimates are handled near the edges. For instance, if we’re at the edge of the search window, do we use zero-padding? Could you point me to the part of the code where this logic is implemented?
I noticed that the code internally includes a threshold on the SNR, but it’s unclear to me if this option is actively used. Additionally, this threshold is not exposed to users. Could you clarify its purpose and usage?

vbrancat · 2025-01-19T01:30:38Z

contrib/PyCuAmpcor/README.md

+cuDenseOffsets.py --reference $reference --secondary $secondary --ww $ww --wh $wh --sw $sw --sh $sh \
+    --mm $mm --kw $kw --kh $kh --gross $gross --rr $rgshift --aa $azshift --oo $oo \
+    --deramp $deramp --outprefix $outprefix --outsuffix $outsuffix \
+    --gpuid $gpuid  --usemmap $usemmap --mmapsize $mmapsize \
+    --nwac $nwac --nwdc $nwdc --workflow $workflow


Is the flag used to switch between the old PyCuAmpcor implementation (i.e., the one following the Fortran code) and the new v2 implementation exposed to the user? If it is, it would be helpful to add it here

vbrancat · 2025-01-19T01:33:11Z

contrib/PyCuAmpcor/examples/cuDenseOffsets.py

+    parser.add_argument('--startpixeldw', dest='startpixeldw', type=int, default=None,
                        help='Starting Pixel down of the reference image.' +
-                             'Default: %(default)s to be determined by margin and search range.')
+                             'Default: %(default)s to use the first pixel.')


I do not understand this. Can you really start from the first pixel if this parameter has not been assigned? Don't you still need a margin equal to the search window?

This is to fix the location of the reference window, for cases to compare results with different window sizes and search ranges. It could lead to out-of-range issues for windows at the edges. Zeros are padded to keep the program running. Warnings are raised to users to alert them of the inaccuracy of these results.

vbrancat · 2025-01-19T01:37:21Z

contrib/PyCuAmpcor/examples/cuDenseOffsets.py

+    parser.add_argument('--wf', '--workflow', dest='workflow', type=int, default=0,
+                        help='workflow (0 = ROIPAC, 1 = GrIMP) (default: %(default)s).')


Could we include this flag in the example? The original Ampcor followed the GRIMP approach but later evolved into the ROIPAC algorithm for computational efficiency. To ensure transparency for users unfamiliar with GRIMP and ROIPAC, we should reframe this in terms of enabling or disabling 2x anti-aliasing oversampling for integer offset detection. Additionally, the help guide should clarify the implications of using or not using 2x oversampling to help users make informed decisions.

vbrancat · 2025-01-19T01:42:47Z

contrib/PyCuAmpcor/examples/cuDenseOffsets.py

+        if inps.use_center != 0:
+            # the given starting coordinate is the center, adjust it
+            startPixelDown -= inps.winhgt//2


My understanding is that PyCuAmpcor uses the top-left corner of a pixel as the reference point (origin) for its pixel coordinates. If this assumption is correct, then to calculate the center of a window relative to its top-left corner, you would need to add half of the window size (in both height and width) to the top-left coordinate.

Explained earlier.

vbrancat · 2025-01-19T01:43:16Z

contrib/PyCuAmpcor/examples/cuDenseOffsets.py

+        # use margin + halfSearchHeight instead
+        startPixelDown = margin + inps.srchgt


This sort of contradicts what you have written in the help i.e., using the first pixel

explained earlier.

vbrancat · 2025-01-19T01:55:53Z

contrib/PyCuAmpcor/src/cuAmpcorParameter.h

+    // correlation surface
+    int2 corrWindowSize; // 2*halfSearchRange + 1
+    int2 corrZoomInSize; // zoomWindowSize+1
+    int2 corrZoomInOversampledSize; // corrZoomInSize * oversamplingFactor


Nit-picking. Since we have multiple stages of oversampling, can we specify which oversampling factor are we referring to here? I assume this refers to the cross-correlation surface oversampling factor?

vbrancat · 2025-01-19T01:56:38Z

contrib/PyCuAmpcor/src/cuAmpcorParameter.h

+    int2 corrZoomInSize; // zoomWindowSize+1
+    int2 corrZoomInOversampledSize; // corrZoomInSize * oversamplingFactor
+
+    int corrStatWindowSize;     ///< correlation surface size used to estimate snr


Suggested change

int corrStatWindowSize; ///< correlation surface size used to estimate snr

int corrStatWindowSize; ///< correlation surface size used to estimate correlation surface statistics (e.g., peak SNR, covariance matrix)

covariance only needs a few pixels around the peak. The statement was correct.

vbrancat · 2025-01-19T02:00:51Z

contrib/PyCuAmpcor/src/cuAmpcorParameter.h

+
+    int corrStatWindowSize;     ///< correlation surface size used to estimate snr

    float thresholdSNR;      ///< Threshold of Signal noise ratio to remove noisy data


My understanding is that PyCuAmpcor computes each offset independently of its associated statistics. If this is correct, why is there a threshold on the SNR? Although this threshold is not exposed to users, what would happen if it were? Does PyCuAmpcor skip the computation of the sub-pixel offset for estimates with an SNR below the threshold? SNR is an unreliable metric for assessing the quality of offset estimates. If outlier removal is necessary, I recommend using the covariance matrix as a more robust alternative.

Deleted. PyCuAmpcor leaves filtering, cullings for users.

vbrancat · 2025-01-19T02:02:04Z

contrib/PyCuAmpcor/src/cuAmpcorParameter.h

    std::string offsetImageName;       ///< Offset fields output filename
    std::string snrImageName;          ///< Output SNR filename
    std::string covImageName;          ///< Output variance filename
+    std::string peakValueImageName;      ///< Output correlation surface peak value filename


Suggested change

std::string peakValueImageName; ///< Output correlation surface peak value filename

std::string peakValueImageName; ///< Output normalized correlation surface peak value filename

Let's specify that is normalized, otherwise the use of this metric is meaningless

vbrancat · 2025-01-21T17:10:43Z

contrib/PyCuAmpcor/src/cuEstimateStats.cu

-            covValue[idxImage] = make_float3(99.0, 99.0, 0.0);
+            covValue[idxImage] = make_real3(99.0, 99.0, 0.0);


To my understanding, if the Gaussian curvature is too small (less than 0.01) then we assign a noData value of 99. Why we do this and where is the threshold of 0.01 come from?

This follows the ROIPAC. We could make it a user configurable parameter.

hfattahi · 2025-01-21T19:00:54Z

@lijun99 Thank you for adding the PR description.
I agree with @vbrancat 's points about the GRIMP/ROIPAC terminology. Perhaps the correct thing is to be just clear about the changes made to the PyCuAmpcor. Main reason is that the base of GRIMP was also original JPL's roipac ampcor. The fact that upfront upsampling has been dropped at some point in some version of it for computational efficiency while kept in other versions such as the one in GRIMP may be just a source of confusion.

We see that one big change contributing to the big PR is the functionality to build the code with both single/double precision. Since our tests have not demonstrated a meaningful difference for the offsets results with the two builds, we are a bit hesitant if we should go that way. Or at least we are not sure if the implementation approach for adding double precision is the most efficient way. Are you open to:

1- remove double precision build from this PR. We expect that would significantly reduce the number of lines of codes in this PR
2- If you definitely want to add double precision and have a good reason for it, let's discuss if templating the code would be a better approach and have a separate PR in future.

vbrancat · 2025-01-23T07:12:01Z

contrib/PyCuAmpcor/examples/cuDenseOffsets.py

  # offset and its geometry
  # tip: re-run with --full/out-geom and without --redo to generate geometry only
-  cuDenseOffsets.py -r ./SLC/20151120/20151120.slc.full -s ./SLC/20151214/20151214.slc.full --outprefix ./offsets/20151120_20151214/offset --ww 256 --wh 256 --sw 8 --sh 8 --oo 32 --kw 300 --kh 100 --nwac 100 --nwdc 1 --gpuid 2 --full-geom ./geom_reference --out-geom ./offset/geom_reference
+  cuDenseOffsets.py -r ./SLC/20151120/20151120.slc.full -s ./SLC/20151214/20151214.slc.full --outprefix ./offsets/20151120_20151214/offset --ww 256 --wh 256 --sw 8 --sh 8 --oo 32 --kw 300 --kh 100 --nwac 100 --nwdc 1 --gpuid 0 --full-geom ./geom_reference --out-geom ./offset/geom_reference


I am not able to comment to the exact part of the code but the current help does not contain all the deramping options i.e.,

0: no-deramping, use magnitude for oversampling
1: oversample complex patches and deramp with linear ramp
2: no-deramping, use complex data for oversampling

Can you update the command line function and the help accordingly? Thanks

…umber of NCC passes

…Offsets.py

…f the oversampled correlation surface

…d to be a factor of #SMs

…ssue

…the warnings generated by -Wreorder

+ self.workflow either 0 (old two-step) or 1 (new one-step) + update this option for both cpu and gpu function + add workflow keyword into applications.topsApp.TopsInSAR

…for correlation peak values

…flow

lijun99 · 2025-06-23T18:19:46Z

I have merged all the pycuampcor updates to a separate repo, which can be embedded as a submodule in isce2/isce3: here is an example. @rtburns-jpl , please let me know whether this is what you prefer. ( I have renamed the module in lower cases, per the conda requirement).

Lijun Zhu added 25 commits May 28, 2022 16:17

PyCuAmpcor: merge isce-framework#622 by rtburns-jpl

cf86b83

pycuampcor: deallocate memory objects

3a4fa40

Merge branch 'main' into pycuampcor

be9d498

pycuampcor: add debug messages for loading the images

1be8be4

PyCuAmpcor V2: start with oversampled images

9ecd6b9

PyCuAmpcor: fix the formula for determining the final offset

e85f962

PyCuAmpcor: oversampling procedure adjustment

5f48f3e

PyCuAmpcor: add extra pads to correlation surface

6704ee9

PyCuAmpcor: use split spectrum for even sequence oversampling

85c8ec1

PyCuAmpcor: add margin to the loaded secondary chip as well

fe15e7a

PyCuAmpcor: add a python script to plot offset fields

3d7de1b

PyCuAmpcor: use a larger area to oversample the reference image

9ad2269

PyCuAmpcor: first attemp with double precision - use CUAMPCOR_DOUBLE …

68cf742

…flag

PyCuAmpcor-DB: first running version

90fe814

PyCuAmpcor-DB: fixed some errors to produce the correct results

bcd2b91

PyCuAmpcor-DP: disable time domain algorithm for (os) window larger t…

0cbd93b

…han 640 with DP

PyCuAmpcor-DP: change default to single precision

67250aa

PyCuAmpcor-v2: solution to time domain code for gpus with 48K shared …

fd28b41

…memory

PyCuAmpcor-V2: in zero padding for FFT oversampling, resort to conven…

c656d72

…tional approach to side Nyquist frequency (N/2) with negative frequency

PyCuAmpcor-v2: add correlation surface statistics (snr and variance) …

1da8c64

…back

PyCuAmpcor-v2: add the correlation surface peak values to output

baf1df3

PyCuAmpcor-v2: rename cuAmpcorChunk to cuAmpcorProcessor and use the …

cc25c3b

…virtual class for different workflows

PyCuAmpcor-v2: bring back ROIPAC/ampcor workflow (set as default), br…

03e3300

…ing back Sinc Oversampler

PyCuAmpcor-v2: repair Readme

a6764b8

hfattahi requested review from gmgunter, rtburns-jpl, vbrancat and xhuang-jpl January 17, 2025 05:37

vbrancat reviewed Jan 21, 2025

View reviewed changes

vbrancat reviewed Jan 23, 2025

View reviewed changes

Lijun Zhu and others added 16 commits January 28, 2025 15:51

PyAcuAmpcor-v2: add debug info to cuderamp=1

7d3eb2c

PyAcuAmpcor-v2: fix the search range in presence of extra pads

a3eb6a4

PyAcuAmpcor-v2: rename two workflows to TwoPass and OnePass wrt the n…

dd21e68

…umber of NCC passes

PyCuAmpcor-v2: clarify notations and notes, per PR comments

5d070d5

PyCuAmpcor-v2: add description for the offset visualization tool plot…

538704f

…Offsets.py

PyCuAmpcor-v2: revert to custom mmap image loader

fa3231e

PyCuAmpcor-v2: 1) fix data types for DP; 2) reduce the search range o…

936d7bc

…f the oversampled correlation surface

pycuampcor-v2: fix an error in determining the corrwindow size

fbcc209

pycuampcor-v2: fix also the corrwindow size error for two-pass workflow

80a7c45

pycuampcor-v2: add get_sm_count method - the batch size is recommende…

cecdee8

…d to be a factor of #SMs

patch getLineBand in some fortran code

68a4199

Merge branch 'main' into pycuampcor-v2-devel

1670d80

Merge branch 'fortran-patch' into pycuampcor-v2-devel

1ea8c27

PyCuAmpcor: patch the scons script for the two different vironments i…

9df1506

…ssue

pycuampcor-v2: change the order of parameters in classes to suppress …

54070c2

…the warnings generated by -Wreorder

pycuAmpcor add new workflow to TopsProc, topsApp

7acaf20

+ self.workflow either 0 (old two-step) or 1 (new one-step) + update this option for both cpu and gpu function + add workflow keyword into applications.topsApp.TopsInSAR

yuankailiu force-pushed the pycuampcor-v2-devel branch from 4638b1a to 7acaf20 Compare April 10, 2025 23:12

Lijun Zhu and others added 5 commits April 10, 2025 20:31

pycuampcor for topsApp: change the slc access method; add the output …

3269450

…for correlation peak values

add OFFSET_PEAKVALUE_FILE to TopsProc

b99a56c

pycuampcor in topsApp: fix the runDenseOffsets script

b9572cd

revert typo changes

d2b4104

pycuampcor-v2: fix the distance in variance estimate in one-pass work…

c417f5f

…flow

hfattahi mentioned this pull request Jul 29, 2025

Pycuampcor V2 isce-framework/isce3#76

Open

		parser.add_argument('--wf', '--workflow', dest='workflow', type=int, default=0,
		help='workflow (0 = ROIPAC, 1 = GrIMP) (default: %(default)s).')

		# use margin + halfSearchHeight instead
		startPixelDown = margin + inps.srchgt

	int corrStatWindowSize; ///< correlation surface size used to estimate snr
	int corrStatWindowSize; ///< correlation surface size used to estimate correlation surface statistics (e.g., peak SNR, covariance matrix)


		int corrStatWindowSize; ///< correlation surface size used to estimate snr

		float thresholdSNR; ///< Threshold of Signal noise ratio to remove noisy data

	std::string peakValueImageName; ///< Output correlation surface peak value filename
	std::string peakValueImageName; ///< Output normalized correlation surface peak value filename

		covValue[idxImage] = make_float3(99.0, 99.0, 0.0);
		covValue[idxImage] = make_real3(99.0, 99.0, 0.0);

Pycuampcor v2 #922

Are you sure you want to change the base?

Pycuampcor v2 #922

Uh oh!

Conversation

lijun99 commented Jan 16, 2025 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

hfattahi commented Jan 17, 2025

Uh oh!

vbrancat left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hfattahi commented Jan 21, 2025

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

lijun99 commented Jun 23, 2025

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

lijun99 commented Jan 16, 2025 •

edited

Loading