Optimizing and Minimizing


In the last post I wrote about my attempt to estimate crater radius and center location using Sartopo data. Lets review.

Screenshot at Jul 25 11-31-12

I use this hypothetical circle to idealize what what a crater would look like. In reality, it isn’t a perfect circle. We can see in topography data that the rims are not equal distances from the center. The point closest to the center is dc, which separates the left and right into two right triangles. Each side can be thought of as some fraction x or y or the total profile length Dp. We can define Dp as follows.

Screenshot at Jul 25 11-35-46

Script based minimizing

I approached this thinking, how can I minimize these two equations? I thought, x, r_left, r_right, and d_c are all solvable if we assume we know where the center is. So I created a script to solve for Dp at each possible point of the center. It seemed promising, where each equation seemed to produce a circle of points that approached the desired Dp/ Subtracting the two from one another, which should produce a 0, would hone in on exactly where they meet, focusing on two points, one near the center found using the radar image and one far outside of view. Unfortunately, it was too far off. I thought it was simply a mistake somewhere in the script, so I proceeded to create a write up starting by reworking through the math I used.

We already have the equation for Dp above, but we also can find the other variables,Screenshot at Jul 25 11-54-20.png

such that Dp=Dp. The calculations will be slightly rounded, but the net result approaches Dp because I use circular logic. This was not what I was getting. I was perplexed, and after staring for a while, Kevin lent me a hand. He pointed out that I was multiplying by x instead of dividing. That concluded this failed attempt.

Minimizing with MATLAB function

Beginning again with Dp,

Screenshot at Jul 25 12-01-32

This still leaves us with an function that doesn’t use Dp. I made a mistake in my logic, which I only just realized now, even after writing it up.

Screenshot at Jul 25 12-08-53

I moved x and y to the other side of the equation, but also kept it on the left side. This lead to incorrect and inconsistent results.

Luckily, the steps I used, of which I have more to say (later in this post) are already in place, so I can easily adjust the function and run it again. Where the new function is found to be as follows.

Screenshot at Jul 25 12-18-01

Keep in mind, x can be solved as a function of the other 3 variables are a matter of geometry. The other three numbers have a more direct relation to the center as a function of Dp. To my surprise, the initial results are more promising.Screenshot at Jul 25 12-27-53.png

Literature put the center at the red dot in the center, and this estimate puts it at the yellow dot with a red rim. There is a bit more I’d like to do such as play with the constraints and plot a circle to see how well it fits the rim and radar image.

Visualizing and translating the results

This is a necessary thing to discuss because it took several days to take these variables and output values for the center lat and long.  I ended up working this out multiple times from different points of view to find the center.

Note: phi is found using the dx and dy between the left and right rim.

I thought the logic was sound, but I kept getting results that did not make sense. I have x and y values for each point on the profile, but when I found the distance between points using the x and y values, they did not translate into lat and long conversions appropriately. The biggest problem is that the y points don’t change as much as the longitude. X appears to change consistently with longitude, but I found it still to be inaccurate. Screenshot at Jul 25 14-17-56

I compared changes in longitude with changes in x. I converted each of these to the other and found the two did not translate. The obvious question to ask is how I am converting between the two. Distances in kilometers were converted to lat and longs with,

rLon = dx/ (R *pi()/180 * cos(deg2rad(latitude)));

%where the inputed latitude for the center  is used

rLat = dy / (pi()/180 * R);

I wondered if there was something about this equation that wasn’t consistent or that didn’t apply for Titan, so I tried to compare it with another method that goes from changes in lat/lon to changes in kilometers. The method comes from one Andrew Hedges. I chose this because its the recommended method from JPL.

dlon = lon2 – lon1
dlat = lat2 – lat1
a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
c = 2 * atan2( sqrt(a), sqrt(1-a) )
d = R * c

I wondered if I had inputed this wrong, so I changed R to that of earth and tested it with the calculator at the site I link to above and it matched. I used this to calculate the profile length and it was off ~11km from the dbidr. Using just along x, I got a number ~5 times higher than the longitude and nearly as much higher than the dbidr value, meaning its an order of magnitude larger. dbidr ~=77km, lat/long~=66km, dx ~=300km. Therefore, lat long more reliable than the x and y values.

I don’t know if I am wrong to use these lat and lon conversions, but the x and y numbers do not work for solving in the x and y domain. So I used the initial lat and long values to to define the left rim, the dc intersection, and the center  in x and y. Changes in lat and long between the left rim and dc, plus the results of the minimization give new x and y for the center. I can then find the change in x and y between dc and the center, or from lr and the center. I did both and got the same resulting lat and lon.

I feel like there must be something about the nature of the sartopo x and y data that I do not understand because the alternative is that it is wrong. Obviously the error is probably with me, but I have gone through this so many times in some many ways to try and hone in on the problem and its lead me here.

Its also worth noting, I wrote the script to use the dx/dy, dlat/dlon depending on the slope of the profile and its general location relative to the center because that will influence whether we are moving to higher or lower lat/lons. Meaning, it isn’t about how I am using the changes in x/y.




July Update

Crater Topography (depths and roughness)

I read through Catherine’s paper on Titan crater topography, and realized that roughness wasn’t the way to go about it. I’ve since realized she left a comment saying as much on my last post. I’ve changed the roughness parameter to relative depth, the measurement that Catherine uses to see how much a crater has changed from its state (as compared to Ganymede craters).

I calculate depths using the same technique that Catherine does finding the lowest point on each side of the center of the crater. I used the depths for Ganymede (Bray and Schenck papers) I found a while back when I worked to create meshes for the tekton code. I’m still faced with the same issues in regards to the very large craters, but that isn’t the main priority at the moment. I know Catherine included Menrva in her data, so I am curious how she determined what a Ganymede depth for a 400+km crater should be.

Best Fit Circle

The circle fit that I discussed before is a function developed by Izhak Bucher in 1991. It works by using a little matrix algebra that I am having a little trouble follow. I had hoped to decipher exactly how it works before I updated you all but wasn’t successful. We start with a matrix A consisting of three columns: x, y, ones(length of x). This is divided into a vector = -x.^2+-y.^2. The result is a vector of the x and y center points and can calculate. My knowledge of matrix division is a little rusty, and from what I could find a\b = inv(a)*b but inv only works with square matrices. Without really understanding this step I can’t really explain the method exactly, merely that it tries to best fit the circle to all the given points.

I just reread it again, and I think I may have misunderstood what you were getting at. It doesn’t use RADAR imagery if that is what you are worried about. It uses the position of the rim from multiple flyovers.

If we look back at one of my figures from my last post, we can see how the rim is identified in each sartopo strip. It’s these sartopo points that the circle is being fitted to. The radar image is merely there as reference.


I thought this was a rather good way to incorporate all the data when we have multiple points available. It allows us to fit the shape to all the data and back out exactly where the center should be. Its when we have less than two fly bys, like in the example I gave last post, where it gets difficult to solve. I don’t like having to make the assumptions, but I’m finding it difficult to find a better way when we have such little data. A circle requires 3 or more points.

Single Sartopo Profile Approach

The use of the Chi-squared test is an interesting approach I didn’t think of. I haven’t done it yet because I hadn’t realized I had a response on my blog. However, I have used it before in a statistical analysis class back at Tech. I still have all my notes for this class, so I think I can do this. However, from what I can ascertain from an initial overview, it compares the modeled to the observed. I assume I would input a range of possible values for each variable to get different fits. Then the chi-squared test would be done on the fit and the actual data. It seems that I am to just compare two numbers.

r_left = sqrt(d_c^2 + (D_p/x)^2)
r_right = sqrt(d_c^2 + (D_p/y)^2)

The variables x and y are not exactly two unless the two radii are the same. So, we have two equations with five unknowns: r_l, r_r, d_c, x, and y. You could try a minimization process to get the ‘best guess’ for these values. Start with reasonable assumptions and see which combination produces the smallest chi-squared value when compared to the known value of D_p.

D_p = x * sqrt(r_left^2 – d_c^2) = y * sqrt(r_right^2 – d_c^2)


Here I am just calculating what the profile length (two calculated values) would be and comparing it to the single value it actually is. Why would I try to do a chi analysis when I could just find the closest value? I don’t understand what I would be performing it on.

Error and Uncertainty

I still have a ways to go with uncertainty. I updated the code to pull the random and systematic errors, but I don’t entirely understand how they fit together. Does one incorporate the other? I’ve added the two together for now and can add them to the outputted result of the code.


Thanks for reading, comments appreciated!