Large-scale nonlinear optimization on a triangular mesh

[Eric Verner]

Hello,

I am trying to optimize a surface represented by a 2D grid of points to produce normal vectors of the surface that align with provided target normal vectors. The grid size is likely to be between 201x201 and 1001x1001. That means that the number of variables will be 40,000 to 1,000,000, as I am only modifying the z-coordinates of the mesh points. I originally posted this on Stack Overflow and was grateful that Dr. Agarwal responded.

Here is the original Stack Overflow post.

http://stackoverflow.com/questions/26336036/how-to-solve-large-scale-nonlinear-optimization-problems-with-ceres

I have been trying to implement a solution suggested by Dr. Agarwal but am running into problems. In my first attempt, I created a Dynamic Autodiff cost function for every mesh point and passed in the adjacent points as constants. This resulted in all of the interior (non-edge and non-corner) points staying at zero because the vertex normals are the average of the face normals, i.e., a change in the height of a vertex would have no change on the normal because all of the face normals would cancel each other out in the average.

The next solution we tried was to pass in the heights of the adjacent points as parameters. The problem with this solution is that the heights within a cost function are optimized, but the parameters in each cost function are not communicating with each other. In other words, each cost function has its own version of what its neighborhood of points should look like, but they all conflict. We collected all of the optimized residuals from each cost function and averaged them, but the solution was not usable.

We would like to have cost functions at adjacent vertices operate on the same data and work together to find a solution for the whole mesh. How do we do this? We have a Mesh object that holds the vertices in a private variable X. How do we pass this data into the cost functions so that the solver finds a global solution and not just individual solutions for each cost function?

I would appreciate help with this, as we are still working on the problem after taking a break for about 2 months.

Thanks,

Eric

[Keir Mierle]

Hi Eric,

The problem you describe should be solvable with Ceres. Do you have recent code to show what you are attempting? I am not sure what you mean by this statement: "the parameters in each cost function are not communicating with each other". Generally Ceres handles this case properly. My interpretation of the problem from the SO post is as follows:

Minimize the difference between a given set of surface normals and the normals from a mesh

Optimizing over the Z-coordinate of the mesh

The mesh starts as a flat patch

This raises a key question: How do you determine the surface normal of a mesh? There are many ways to do it- you could take a 1-neighborhood around each normal sample and fit a plane to the z-heights, then take the normal to the plane. You could do a sophisticated quadratic fit to any number of surrounding points, then take the tangent of the quadratic.

With that said, here is how I would try to solve the problem:

For each normal sample where we have a parameter (assume interior for now)

  Create a cost function taking 9 parameters - the current sample & neighbors

    For each neighbor iterating *around* the center, compute a surface normal formed by the neighbor, the current sample, and the previous neighbor - store it

    Average the 9 surface normals to get a normal estimate for the current point

    Subtract the computed normal from the expected normal to form a 3-dimensional residual

This means you will have to create 5 cost functions implementations - one for each top/bottom/left/right/interior - due to the structure of Ceres' autodiff cost function. Hopefully you can reuse the normal computation code by factoring it out.

Alternately, you could do something sneaky like create a 1-width border around your mesh, then set the border parameters to constant with Ceres. You could then use the interior-only cost function everywhere. That would save you the hassle of making 5 implementations at a minor increase to computation cost (since Ceres is smart about handling fixed parameters).

Happy optimizing,

Keir

--
You received this message because you are subscribed to the Google Groups "Ceres Solver" group.
To unsubscribe from this group and stop receiving emails from it, send an email to ceres-solver...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/e57da29b-e68f-42b9-9b45-4137324d129b%40googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Keir Mierle]

On Tue, Dec 16, 2014 at 9:17 PM, Keir Mierle <mie...@gmail.com> wrote:

Hi Eric,

The problem you describe should be solvable with Ceres. Do you have recent code to show what you are attempting? I am not sure what you mean by this statement: "the parameters in each cost function are not communicating with each other". Generally Ceres handles this case properly. My interpretation of the problem from the SO post is as follows:

Minimize the difference between a given set of surface normals and the normals from a mesh

Optimizing over the Z-coordinate of the mesh

The mesh starts as a flat patch

This raises a key question: How do you determine the surface normal of a mesh? There are many ways to do it- you could take a 1-neighborhood around each normal sample and fit a plane to the z-heights, then take the normal to the plane. You could do a sophisticated quadratic fit to any number of surrounding points, then take the tangent of the quadratic.

With that said, here is how I would try to solve the problem:

For each normal sample where we have a parameter (assume interior for now)

  Create a cost function taking 9 parameters - the current sample & neighbors

    For each neighbor iterating *around* the center, compute a surface normal formed by the neighbor, the current sample, and the previous neighbor - store it

    Average the 9 surface normals to get a normal estimate for the current point

    Subtract the computed normal from the expected normal to form a 3-dimensional residual

This means you will have to create 5 cost functions implementations - one for each top/bottom/left/right/interior - due to the structure of Ceres' autodiff cost function. Hopefully you can reuse the normal computation code by factoring it out.

Alternately, you could do something sneaky like create a 1-width border around your mesh, then set the border parameters to constant with Ceres. You could then use the interior-only cost function everywhere. That would save you the hassle of making 5 implementations at a minor increase to computation cost (since Ceres is smart about handling fixed parameters).

Upon further thought, the border idea is a bad since the edges would get pinned down. You'll have to create the 5 separate cost functions, unfortunately.

[Eric Verner]

Hi Keir,

The way we find the surface normals is pretty standard. For every point, we compute the normal vector of each adjacent face, average them, and normalize them. The adjacent points are provided by a Mesh object, which is passed into the function. We only need one cost function to implement this. The procedure goes like this:

1) Compute normal vector of each adjacent face.

2) Normalize the face normal vectors to have magnitude of 1.

3) Average the face normal vectors to find the vertex normal vector.

4) Normalize the vertex normal vector to have magnitude of 1.

Here is the objective function, which performs this calculation:  http://pastebin.com/bXkgkzNB

Here is the code in the main function that solves the optimization problem using Ceres: http://pastebin.com/AxUdVCA6

Here is the member function in the Mesh object that provides the adjacent points to the objective function: http://pastebin.com/BgyDX8qa

The problem is that the multiple residual blocks don't seem to work together to produce a good solution.

Thanks for your help.

Eric

[William Rucklidge]

The code you posted isn't going to work very well. Your cost function is computing the normal difference based on a single z value from Ceres, and the existing z values in the Mesh object. But Ceres is trying to optimise all the z values simultaneously. As it stands, your function can compute \sum_i normaldiff_i(z_i) where normaldiff_i uses (other z_i) than the values in the existing Mesh... but there are two problems

1. The values in the existing Mesh are not the values that Ceres is currently trying out in its search of parameter space

2. Your function can only return d normaldiff_i / dz_i - the scalar derivative of the residual for point i wrt z_i. But normaldiff_i also depends on z_j for j \in N(i) - all of i's neighbours... so the Jacobian that Ceres computes will be incorrect.

What you want to do is something more like (this is a more explicit version of what Keir suggested, but I'm writing it out only for the non-edge, non-corner case):

void operator()(double* z_me, double* z_nw, double* z_n, double* z_ne, double* z_e, double* z_se, double* z_s, double* z_sw, double* z_w, double* residual) {

   ... compute residual as a function of these 9 z values ...

}

Then you'd construct your cost functions something like (still ignoring the edges):

typedef ceres::AutoDiffCostFunction<MyCostFunctor, 3, 1, 1, 1, 1, 1, 1, 1, 1, 1> MyCostFunction;

    for (int m = 0; m < NV; m++ ) {

      double* z_me = ... get mutable pointer to appropriate z value ...;

      double* z_nw = ... and likewise for its NW neighbour, and furthermore for all the other neighbours;

      problem.AddResidualBlock(new MyCostFunction(new MyCostFunctor), NULL, {z_me, z_nw, z_n, ...});

}

You could use a DynamicAutoDiffCostFunction and pass in different sizes for the different cases (edge, corner, middle) - but then your cost function has to know what case it is, and has to use the z values passed into it, not any living in the Mesh.

A couple of other notes: in auto-diff functions, you need to be careful when you use sqrt() as its derivative has a singularity at zero. So if your normal vectors will ever be zero (e.g., if your mesh starts out as flat), then you'll have trouble; best to handle that as a special case. Also, you don't need to compute the sum of squares of the three components of your residual vector - just put the three components into the three elements of *residual (Ceres knows there should be 3 because of the "3" in the templating of AutoDiffCostFunction) and it will compute the squared magnitude for you. As your code is written, you're returning the squared magnitude, which Ceres will then square again, which is not what you want.

-wjr

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/4a69f88e-5d3f-483c-aad0-d8ab1edebd75%40googlegroups.com.

[Eric Verner]

Thank you, wjr, for your help! I will try your recommendations out.

[Francesco Callari]

[ Coming late to this party, with a general suggestion from experience... ]

The approach you indicate has some problems - if there is nothing constraining the geometry of the solution other than the prescribed normals, you can easily get local minima with unreasonable topologies - e.g. facets shrinking to a point or a sliver, or overlapping each other.

 

The problem you describe is a rather common one in computational fluid dynamics. If the mesh in question is a boundary surface, the usual approach is to embed it into a level set and evolve it with the given normals as a boundary condition. The "evolution" is defined by a process appropriate for the case at hand (advection, convection, inflation, ...) and is constrained by appropriate physical properties (e.g. body elasticity).

If your initial configuration is expected to be "very close" to the optimal solution, for example because the optimal solution is just a "smoothed" version of the initial state, you can try hacking it. An approach I have used for meshes representing human faces is a simple variant of gradient descent, where the "gradient" is defined at each facet by the component of (prescribed_normal - current_normal) that is orthogonal to current_normal. Each of these components is then interpreted as producing a torque at the center of the facet, and the vertices are then moved by a small multiple of the sum of the motions implied by the torques. Iterate and pray...

HTH

Franco

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/0a90c5ab-6d60-46d6-8307-ee34c1f3b251%40googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--

Franco Callari <fgca...@gmail.com>

            EC67 BEBE 62AC 8415 7591  2B12 A6CD D5EE D8CB D0ED

I am not bound to win, but I am bound to be true. I am not bound to succeed, but I am bound to live by the light that I have. (Abraham Lincoln)

[Yuliy Schwartzburg]

Hi everyone,

I have quite a lot of experience with this type of problem. Franco, as Eric is only optimizing the z values of a heightfield, the facets will not shrink or overlap. Hacking it is not necessary.

Eric, wjr is right, right now you pass the entire mesh to your cost function and look up the adjacency every time, but your connectivity is fixed for the optimization (and it should be anyway for tractability). You can indeed use a DynamicAutoDiffCostFunction to handle all the cases together, and send in all the parameters in the one ring neighborhood of the point when you build the cost function at the start.

Best,

Yuliy

[Eric Verner]

Hi,

Thanks for your advice, everyone! Real quick, how do I do this?

double* z_me = ... get mutable pointer to appropriate z value ...;

Would something like this work?

double* temp = new double[3];
temp = mesh.vertex(m);
double* z_me = temp;

Where inside the Mesh class...

double** X_;     // Vertex coordinates [NV x 3]
...
double* vertex(int m) { return X_[m]; }

Additionally, is it possible to just create an array of mutable z coordinates (e.g., Z), which would allow me to replace (1) with (2)?

(1) problem.AddResidualBlock(new MyCostFunction(new MyCostFunctor), NULL, {z_me, z_nw, z_n, ...});

(2) problem.AddResidualBlock(new MyCostFunction(new MyCostFunctor), NULL, Z);

This would save time in writing extra cases.

Thanks,
Eric

[William Rucklidge]

On Sun, Dec 21, 2014 at 6:13 PM, Eric Verner <eve...@gmail.com> wrote:

Hi,

Thanks for your advice, everyone! Real quick, how do I do this?

double* z_me = ... get mutable pointer to appropriate z value ...;

Would something like this work?

double* temp = new double[3];
temp = mesh.vertex(m);
double* z_me = temp;

That would leak the value allocated with "new double[3]".

 

Where inside the Mesh class...

double** X_;     // Vertex coordinates [NV x 3]
...
double* vertex(int m) { return X_[m]; }

Assuming that each vertex coordinate is {x, y, z}, then you can do

double* z_me = mesh.vertex(m) + 2; (or equivalently &(mesh.vertex(m)[2]))

The goal is to get a pointer to the place where you want the final value Ceres computes to be written. It is crucial that different cost functions that use this same value use the same pointer - Ceres needs to know that one vertex's z_me is the same as the z_w for the one to the east, and so on.

Additionally, is it possible to just create an array of mutable z coordinates (e.g., Z), which would allow me to replace (1) with (2)?

(1) problem.AddResidualBlock(new MyCostFunction(new MyCostFunctor), NULL, {z_me, z_nw, z_n, ...});

(2) problem.AddResidualBlock(new MyCostFunction(new MyCostFunctor), NULL, Z);

Do you mean something like having a single Z array containing (usually) 9 values? No, that's not possible - basically you've shattered your grid into a bunch of independent points. How can Ceres know that cost function A's third Z value is the same value (height as the same grid point) as cost function B's seventh?

You can use one of the "Dynamic" flavours of cost functions so that sometimes you pass in a 9-vector of Z pointers, sometimes a 6-vector, sometimes a 4-vector. But your cost function internally needs to know whether it's an interior, edge, or corner point (and which edge/corner it's on) so that it can properly use the values Ceres hands it.

-wjr

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/12d0fa9f-f126-45c0-b9c5-da0c9d2809cb%40googlegroups.com.

[Eric Verner]

Thanks for the tip on the mutable coords.

"Do you mean something like having a single Z array containing (usually) 9 values? No, that's not possible - basically you've shattered your grid into a bunch of independent points. How can Ceres know that cost function A's third Z value is the same value (height as the same grid point) as cost function B's seventh?"

What if Z is an array of pointers to different points in the mesh? Wouldn't Ceres be able to function correctly as long as each cost function was given the correct memory addresses to the data? It is not necessary to separate them into different variables in the cost function call and declaration, right?

"You can use one of the "Dynamic" flavours of cost functions so that sometimes you pass in a 9-vector of Z pointers, sometimes a 6-vector, sometimes a 4-vector. But your cost function internally needs to know whether it's an interior, edge, or corner point (and which edge/corner it's on) so that it can properly use the values Ceres hands it."

Given that the connectivity pattern is uniform throughout the mesh, can we give the Dynamic cost function pointers to the mesh points in the same order every time (e.g., CW or CCW), and then just tell the cost function to iterate in one direction through the array? That seems like it would work without having to specify the different cases (4 corners, 4 edges, 1 interior).

BTW, thanks for all your help, wjr. I really appreciate it.

-Eric

[Eric Verner]

Actually, with a Dynamic cost function, each parameter would be added separately, not with an array. The inputed array of pointers would only be necessary for a non-Dynamic cost function. However, the Dynamic cost function would see an array, not separate variables.

[William Rucklidge]

On Mon, Dec 22, 2014 at 11:37 AM, Eric Verner <eve...@gmail.com> wrote:

Thanks for the tip on the mutable coords.

"Do you mean something like having a single Z array containing (usually) 9 values? No, that's not possible - basically you've shattered your grid into a bunch of independent points. How can Ceres know that cost function A's third Z value is the same value (height as the same grid point) as cost function B's seventh?"

What if Z is an array of pointers to different points in the mesh? Wouldn't Ceres be able to function correctly as long as each cost function was given the correct memory addresses to the data? It is not necessary to separate them into different variables in the cost function call and declaration, right?

An array of pointers would work in concept, and in fact is how the raw CostFunction interface works ("double const* const* parameters" - it passes you an array of pointers). But if you're using (say) AutoDiffCostFunction, then you get a one function parameter per parameter block. Sometimes one mode is better-suited to your problem structure, sometimes the other. For example, for problems like bundle adjustment, you might write a cost functor like

void operator()(const T point[3], const T pose[6], const T intrinsics[3], T residual[2])

since the parameter blocks are segmented into different types and you always get one of each type. Your parameter blocks are more uniform, and you sometimes get different numbers of them, so passing arrays of pointers matches the problem structure better - which is exactly what DynamicAutoDiffCostFunction does.

"You can use one of the "Dynamic" flavours of cost functions so that sometimes you pass in a 9-vector of Z pointers, sometimes a 6-vector, sometimes a 4-vector. But your cost function internally needs to know whether it's an interior, edge, or corner point (and which edge/corner it's on) so that it can properly use the values Ceres hands it."

Given that the connectivity pattern is uniform throughout the mesh, can we give the Dynamic cost function pointers to the mesh points in the same order every time (e.g., CW or CCW), and then just tell the cost function to iterate in one direction through the array? That seems like it would work without having to specify the different cases (4 corners, 4 edges, 1 interior).

Depends on how your mesh class is structured, I guess, and it sounds like it could be made to work if your problem doesn't care if it's given 

z_me, z_n, z_ne, z_e

(for the SW corner)

z_me, z_e, z_se, z_s

(for the NW corner)

If it performs the same computation for these cases, then I guess it just needs to know if it's a corner, edge or interior point.

-wjr

 

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/a8a8a5e8-c754-4f87-9aa1-f30976d8d363%40googlegroups.com.

[Eric Verner]

Hi,

We created an implementation with a DynamicAutoDiffCostFunction and a vector of pointers to the data. In the first iteration, the program finds a solution that is pretty close to the correct answer. However, it continues on for 10 more iterations, and each time, the z values that it is optimizing decay until the change in the parameters is less than the stopping criterion. Additionally, at the end, the values that were found in the optimization are not saved in the addresses provided by the pointers. Does anyone know what is going on or how to fix it?

Here is the current main function: http://pastebin.com/heGNtztQ

Thanks a bunch,
Eric

[Eric Verner]

Sorry, this information is probably necessary to understand the implementation:

double** Mesh::getMutableZCoords(int m) {

  DEBUG_PRINT(("getMutableZCoords: %u\n", m));

  double** zcoords = new double*[num_adjacent_pts(m)+1];
  int* indices = adjacent_pt_inds(m);

  // Copy over the z coordinates using the index array for reference
  zcoords[0] = &X_[m][2];
  for (int i=1; i<num_adjacent_pts(m)+1; i++) {
    DEBUG_PRINT(("indices[%u]: %u\n", i-1, indices[i-1]));
    zcoords[i] = &X_[indices[i-1]][2];
  }

  return zcoords;
}

The () operator in the cost function hasn't changed except that there are three residuals instead of one now.

for ( int p = 0; p < 3 ; p++ ) { e[p] = N[p] - Ntgt_[p]; }

[Sameer Agarwal]

Solver log please

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/33c9f447-5403-4faa-af0d-a1eb6542dca7%40googlegroups.com.

[Eric Verner]

iter      cost      cost_change  |gradient|   |step|    tr_ratio  tr_radius  ls_iter  iter_time  total_time
   0  2.650208e-01    0.00e+00    1.96e+00   0.00e+00   0.00e+00  1.00e+04        0    1.42e-03    1.55e-03
   1  2.650208e-01   -9.96e-01    0.00e+00   2.33e-01  -3.90e+00  5.00e+03        1    9.13e-05    1.74e-03
   2  2.650208e-01   -9.96e-01    0.00e+00   2.33e-01  -3.90e+00  1.25e+03        1    2.90e-05    1.81e-03
   3  2.650208e-01   -9.95e-01    0.00e+00   2.33e-01  -3.90e+00  1.56e+02        1    2.60e-05    1.87e-03
   4  2.650208e-01   -9.94e-01    0.00e+00   2.32e-01  -3.89e+00  9.77e+00        1    2.51e-05    1.93e-03
   5  2.650208e-01   -9.77e-01    0.00e+00   2.15e-01  -3.85e+00  3.05e-01        1    2.51e-05    2.00e-03
   6  2.650208e-01   -8.85e-01    0.00e+00   6.46e-02  -6.96e+00  4.77e-03        1    2.45e-05    2.06e-03
   7  2.650208e-01   -8.84e-01    0.00e+00   1.40e-03  -2.72e+02  3.73e-05        1    2.42e-05    2.12e-03
   8  2.650208e-01   -8.85e-01    0.00e+00   1.10e-05  -3.45e+04  1.46e-07        1    2.70e-05    2.23e-03
   9  2.650208e-01   -8.85e-01    0.00e+00   4.29e-08  -8.82e+06  2.84e-10        1    2.49e-05    2.30e-03
  10  2.650208e-01   -8.85e-01    0.00e+00   8.39e-11  -4.52e+09  2.78e-13        1    2.39e-05    2.36e-03
  11  2.650208e-01   -8.85e-01    0.00e+00   8.19e-14  -4.63e+12  1.36e-16        1    2.35e-05    2.42e-03

Solver Summary (v 1.10.0-lapack-suitesparse-openmp)

                                     Original                  Reduced
Parameter blocks                            9                        9
Parameters                                  9                        9
Residual blocks                             9                        9
Residual                                   27                       27

Minimizer                        TRUST_REGION

Sparse linear algebra library    SUITE_SPARSE
Trust region strategy     LEVENBERG_MARQUARDT

                                        Given                     Used
Linear solver          SPARSE_NORMAL_CHOLESKY   SPARSE_NORMAL_CHOLESKY
Threads                                     1                        1
Linear solver threads                       1                        1

Cost:
Initial                          2.650208e-01
Final                            2.650208e-01
Change                           0.000000e+00

Minimizer iterations                       11
Successful steps                            0
Unsuccessful steps                         11

Time (in seconds):
Preprocessor                            0.000

  Residual evaluation                   0.000
  Jacobian evaluation                   0.001
  Linear solver                         0.000
Minimizer                               0.002

Postprocessor                           0.000
Total                                   0.003

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

[Keir Mierle]

Hi Eric,

It appears you have not included the cost function implementation in the paste. I am also curious why you are running the Ceres optimization in a loop.

Also it appears you are leaking memory:

1.  int* api = new int[mesh.num_adjacent_pts(m)];
2.  api = mesh.adjacent_pt_inds(m);

This creates an array of int pointers on line 1, then promptly overwrites the only pointer to it on line 2, leaking the allocated array from line 1.

In this line it's not clear to me why the loop limit goes one past num_adjecent_pts():

    for (int i = 0; i < mesh.num_adjacent_pts(m)+1; i++) {

Can you post a more complete program? Perhaps if you could make something we could execute, it would be easier to fix.

Thanks,
Keir

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/35278cc7-96af-4687-8c0f-8087010871ca%40googlegroups.com.

[Eric Verner]

Hi Keir,

Here is the current cost function: http://pastebin.com/vxjkAQXQ

Here is the header from the Mesh class, which we use to store the data: http://pastebin.com/52rFAvqZ

We are running Ceres in a loop because we thought it might be better able to find the correct answer. It wasn't...

mesh.num_adjacent_pts(m) returns the number of points adjacent to vertex m. mesh.num_adjacent_pts(m)+1 is used because the pointer to the z value at vertex m also needs to be added to the cost function.

I could post other code here, but binary data files are required to run it. If we can't resolve this problem with the information already provided, I'll post the files here.

Thanks very much for your help!

Eric

...

[Keir Mierle]

Hi Eric,

I believe you are indexing into the parameter vector ("z") incorrectly in your cost function. The parameters are organized as first parameter, then dimension of each parameter. For example, if your cost function took 2 parameters: the first had 2 dimensions (x,y), and the second had 3 dimensions (x, y, z), then you could access them as

z[0] <-- Pointer to first parameter (2 dimensions)

z[0][0] == x

z[0][1] == y

z[0][2] == <---- INVALID! First parameter only has 2 dimensions

z[1] <-- Pointer to second parameter (3 dimensions)

z[1][0] == x

z[1][1] == y

z[1][2] == z

In your case, the cost function takes N parameters of size 1, one for each adjacent vertex, so you should be indexing them as z[i][0] or *z[i]. The indexing order matters because the size of a pointer is different than the size of a Jet, which is what T is during autodiff expansion.

Hope that helps,

Keir

To view this discussion on the web visit https://groups.google.com/d/msgid/ceres-solver/7afeb664-bcee-4482-848a-895fd4ffe9fb%40googlegroups.com.

[Eric Verner]

Hi Keir,

Yea, that did the trick. I used *z[k] instead of z[0][k]. Works like a charm now. Thanks so much for your help. Really generous use of your time.

Eric

...