Mastering ACCIV
In the Section on [SettingUpInput Setting up Input], we discussed guidelines for setting the major parameters for controlling ACCIV. However, we will usually need to change parameters several times, after looking at the output of an ACCIV pass, to extract the best possible velocity fields.
There are a couple approaches to this process. Mike has very fast ACCIV execution on his workstation, so he prefer to do a full (several-pass) run of ACCIV, examine parameters, and then make changes before the next run. Organizationally, this is easiest for him, because he can keep all the “run A” work files in a “work_a” directory, and so on, allowing him to keep track of each run and compare results of different runs.
Xylar advocates a different pass-by-pass approach to changing parameters. With this approach, the user always looks carefully at the results of each pass to decide whether the next step is changing parameters and re-trying or running the next pass. This passwise approach is more computationally efficient, since a full run with mistakes in an early pass will waste computer time running the later passes.
In the Section on Setting Up Input, we discussed guidelines for setting the major parameters for controlling ACCIV. However, we will usually need to change parameters several times, after looking at the output of an ACCIV pass, to extract the best possible velocity fields. This section gives a few tips about how to adjust parameters based on examining ACCIV output.
If the correlation box is too large, an overly smooth velocity field will be produced. This is not good because it loses fine structure such as jet peaks. If the correlation box is too small, “bald” patches will be produced where correlations could not be found. Bald patches can be easily seen in quiver plots of the outScatteredVelocity.h5 file. But if you only look in the outGridVelocity.h5 file, you will not see these bald patches because the gridding scheme interpolates over them.
If the search range is too large, false correlations may be found at large distances. Many of these will be identified by ACCIV as outliers, but the scatter in correlation velocity uncertainty vectors will be greater. The histogram of correlation velocity uncertainties (see figure below) can be useful for diagnosing this. If the distribution of uncertainties is very broad, with a long tail extending to high velocities, the problem may be that the search range is not optimized. A nicely-peaked uncertainty distribution, like the plot below, suggests that the search range and other parameters are good.
Histogram of uncertainties in the velocity field, using the correlation velocity uncertainty metric. This is derived from correlationVelocityResiduals and correlationLocationResiduals in the combinedCorrelationTiePoints.h5 file. A nicely peaked distribution like this one suggests that the parameters are well optimized, but if the distribution is not sharply peaked, and has a long tail with high uncertainties, the search range may be too large.
The stride mainly controls the density of the output velocity field. If the stride is too large, the outScatteredVelocity.h5 file will tend show velocity vectors in a grid-like pattern, and `outGridVelocity.h5' will be too smooth. If the stride is too small, it will take a long time to get a velocity that is just as noisy as if you had used a larger stride. So stride should be optimized to make sure you have a nice density of velocity vectors, compared to the spatial scale of noise in the data.
This parameter controls the smoothness of the velocity field on the grid. The smooth velocity fit, based on Lee et al. (1997) (See Optimizing Parameters), uses a hierarchy of B-splines of increasing resolution. At each level of the reconstruction, the velocity reconstruction at a given control point is altered only if sufficiently many scattered velocity vectors are found nearby (within a 4x4 b-spline patch). The parameter smoothFitMinControlPointScatteredNeighbors controls how many points are required: The larger this number is, the more neighbors must be present at a given control point, and therefore the smoother the velocity reconstruction will be. We have found that a value between 16 and 32 produces reasonable results in most cases, especially for later ACCIV passes. Typically, a combination of larger stride and larger smoothFitMinControlPointScatteredNeighbors in earlier ACCIV passes can lead to a smooth approximation to the velocity field that can be refined in later passes by decreasing either or both values.
Many of the minor ACCIV parameters (those not included in the list of 3 above) have been discussed in previous Sections. Here are some additional notes on a few of them:
minimumTimeSeparation, maximumTimeSeparation: These parameters set rules about which image pairs will be created. These are most useful in passes with short time separations (maps taken within a single planetary rotation). You can use these parameters to skip image pairs with very short separations (a few minutes) that will not be able to detect any useful displacements, but allow longer separations (tens of minutes or an hour).
In longer passes, like with these settings,
earlierImageIndices = [1,2,3,4,5]
laterImageIndices = [6,7]
the TimeSeparation parameters are not needed.
pathVectorCount: This controls the amount of vectors that are added along curved streamlines. These are not really independent vectors, but they improve the velocity field. For images separated by long times (and therefore with path lengths between the head and tail of a tie point that are long), this number should be set to something like 12 or 13 to avoid artifacts in grid velocity field. For images separated by short times, this number can be smaller but there is no harm in having more vectors along each path.
removeScatteredOutliers, outlierThresholdConstant: ACCIV removes outliers as part of its internal, iterative process for finding a smooth fit to the velocity data on a regular grid. By default, it removes points that deviate from the previous iteration of the smooth fit by more than outlierThresholdConstant=8 times the median deviation. If outliers are to be removed (removeScatteredOutliers = true), ACCIV first iteratively constructs a smooth velocity field and scattered velocity vectors along curved paths with outliers included, then removes outliers, and then iteratively constructs a new smooth and scattered velocity field. If a user finds that ACCIV is producing too many outliers, the first strategy should be to try to reduce the search range as much as possible without removing real velocities. If this still doesn't do the trick, decreasing outlierThresholdConstant may be worth a try.
The sensitivity of the final output velocity field to the input parameters used is a serious outstanding issue. It can be difficult to tell when you’ve created the best possible velocity field. Xylar and Mike believe that the best way to quantify systematic uncertainties is to extract velocity fields with different parameter sets, and then compare the results. Mike's intuition is that this process will reveal uncertainties considerably larger than those calculated by the correlation velocity uncertainty technique.
The correlation velocity uncertainty technique does not detect systematic uncertainties caused by the choice of input parameters. Also, calculating this uncertainty produces a full velocity-uncertainty field with a range of values that changes with position on the map. Quoting the median of these values does not accurately capture the uncertainty that may be relevant to a particular part of the final velocity field.
Alternative approaches might include producing a full uncertainty map that shows the spatial variation of the uncertainty, or showing histograms of the uncertainties, since their distribution is decidedly non-gaussian and therefore not well characterized by a mean and standard deviation.
During an ACCIV clinic held in 2011, it became apparent that the execution speed of ACCIV, for identical datasets and control parameters, varied considerably among software and hardware platforms. The fastest execution seems to take place on my Mac Pro workstation using the g++ compiler. Since the clinic, the smooth-fit operation has been completely rewritten and code performance appears to have improved dramatically on all systems we have tested. However, performance may still vary across systems.
There had been talk of parallelizing ACCIV at some point in the past. For very large datasets, this may be useful. However, we have found that even very large data sets can be decomposed into smaller, overlapping subimages. Then ACCIV can be performed on each set of subimages independently, and finally the scattered data can be combined trivially into a single data set. (This final data set may still be too large to perform a smooth fit to the end result, but smooth fits of subsets of the data can be combined if this is desired. We took this approach in creating a velocity map over most of Jupiter using Cassini images.)
Techniques for turning planetary images into map projections are beyond the scope of this ACCIV training. However, map images may be supplied to you which contain small navigation errors. The figure demonstrates how these errors may affect the data:
Schematic of an image footprint on the spheroidal coordinate grid for Jupiter. If the assumed navigation (grey box) is offset from the true navigation (dashed black box), then there will be errors in the map data: longitudes will be erroneously small for all points in the map, but the error will be larger on the left side than on the right side.
ACCIV provides the capability to statistically correct for small navigational errors. The principle is simple: a single pass of ACCIV is run on the original images, and we make the assumption that over the entire map field, all the displacements should average to zero. In that case, we can use the average of the offsets in the correlationTiePoints_i_j.h5 files to determine the mean displacement between images. We then determine the shifts for each image that would minimize the average of the displacements found by ACCIV, and apply these shifts to the data to create navigation-corrected images.
Ideally, the shifts should be used to re-navigate the original images and create new maps. In practice, the original navigation is usually sufficiently good that we can just apply translational shifts to the maps.
Xylar has Matlab code and Mike has IDL code to do this. Here is Xylar’s Matlab code with notes:
Run ACCIV on all pairs of images separated by short times (less than one rotation, and with some overlap, so less than about 2 or 3 hours) For each image pair (i,j) where i < j, compute a "pixel offset" and add entries to the n x n matrix A and the n x 2 matrix b as follows:
pixelOffset(1) = median(lons_j-lons_i)/deltaLon
pixelOffset(2) = median(lats_j-lats_i)/deltaLat
A(i,j) = A(i,j) + 1;
A(j,i) = A(j,i) + 1;
A(i,i) = A(i,i) - 1;
A(j,j) = A(j,j) -1;
The ith diagonal should have minus the number of pairs involving image i, and entry i,j should be 1 if there exists such an image pair, 0 if not.
b(i,1) = b(i,1) + pixelOffset(1)
b(i,2) = b(i,2) + pixelOffset(2)
b(j,1) = b(j,1) - pixelOffset(1)
b(j,2) = b(j,2) - pixelOffset(2)
where lons_i, lats_i, lons_j and lats_j are the contents of /x1, /y1, /x2 and /y2, respectively, in the appropriate hdf5 file (correlationTiePoints_i_j.h5) and deltaLon and deltaLat are the size of one pixel in longitude and latitude, respectively.
Here's an alternative way to build the matrix that Mike had better luck with in cases where you have vastly more tie points for some image pairs than for others:
weight = length(lons1);
pixelOffset = [median(lons2-lons1)/deltaLon, median(lats2-lats1)/deltaLat];
A(i,j) = A(i,j) + weight;
A(j,i) = A(j,i) + weight;
A(i,i) = A(i,i) - weight;
A(j,j) = A(j,j) - weight;
b(i,:) = b(i,:) + weight*pixelOffset;
b(j,:) = b(j,:) - weight*pixelOffset;
With this version, each image pair gets a weight proportional to the number of tie points that were found. As before, replace the last row of A with all ones and the last row of b with zeros.
set the average offset to zero, that is, set the last row of A to all ones and both entries in the last row of b to zero:
A(n,:) = 1
b(n,:) = 0
The best fit offset for all images is computed by solving the linear system A*pixelOffsets = b, (i.e., pixleOffsets = A^-1^*b), which is easily done in Matlab, and in IDL.
Go through all images:
interpolate the image and mask using the new, offset indices write the image to a new file