🚀 go-pugleaf

(Copy emailed to poster.)

I am somewhat familiar with some projects trying to optimize data
sets of brain imaging, etc.  However, despite what you likely would
like to hear, the answers to most of your queries are in fact
somewhat general.  This is simply because nonlinear sytems typically
are nontypical, and any sampling technique likely (not always) will
have to be tuned to make it robust and/or efficient for a given
class of data.  My own optimization code, Adaptive Simulated
Annealing (ASA), can be obtained at no charge from my archive in
the header and signature.
          %A L. Ingber
          %T Adaptive Simulated Annealing (ASA)
          %R [http://www.ingber.com/ASA-shar, ASA-shar.Z, ASA.tar.Z,
          ASA.tar.gz, ASA.zip]
          %I Lester Ingber Research
          %C McLean, VA
          %D 1993

Other sources of codes for nonlinear optimization can be obtained from
http://www.skypoint.com/subscribers/ashbury/nonlinear-programming-faq.html
http://plato.la.asu.edu/guide.html
http://solon.cma.univie.ac.at/~neum/glopt.html

In my archive is a file, asa_papers, giving some references to
projects that have used ASA and who have made some comparions to
other optimization codes on "non-toy" problems.  I think the
conclusion is as I have just stated:  On a hard optimization problem,
you very often have to bite the bullet and get pretty familiar with
the optimization code you choose.

In article <32F685BE.2781E494@bmrl.med.uiuc.edu>,
Jeff Tsao  <jtsao@bmrl.med.uiuc.edu> wrote:
:Hi everybody,
:
:	I've been working on coregistering medical images for quite
:some time now, and have never found a satisfactory answer to the
:following problems, so I always resorted to quick hacks. I've also
:noticed that other researchers have created their own hacks as well.
:Now, I'm revisiting the problems. Hopefully with the help of all the
:experts on these Newsgroups (e.g. computer scientist for optimization,
:pattern recognition, etc., game programmer for 3D graphics, etc.),
:perhaps we can find a definitive solution. I appologize for the length
:of this posting, but I want to make the posting clear to people with
:different academic backgrounds.
:
:+---------------------------------------------------------------------+
: NOTE: This is a serious posting for solving a practical and quite
:       general problem in image registration, so please *DO NOT* turn
:       it into another "RE: New 3D Rendering Method" discussion in
:       comp.graphics.algorithm.
:       Also, please be specific in your response (preferably citing
:       references or source codes for your responses).
:       e.g. Don't say: there is no general answer to such and such...
:            Cite an example and explain your solution.
:       Thanks in advance.
:+---------------------------------------------------------------------+

Most of what followed is a description of your "cost function."
Optiimization generally requires the researcher to well model
his/her system.  Generally, this is independent of the optimization
algorithm per se, but sometimes it is not.  For example, in some
cases using some optimization algorithms, you can make the cost
function more efficient to optimize by _not_ using penalty functions
for constraints, but obeying these constraints by rejection sampling
within the optimization process.

The best guide you often can get to make your optimization more
robust and/or efficient is to appreciate that the output of the
optimization process can be viewed as a "gray box" interface with
your system.  That is, look carefully at some sample runs, printing
out a lot of information.  See if and how local trapping might be
occurring; this often is the best guide on how to try to tune your
optimization code of choice.  Unfortunately, there most often will
be little guidance gained, but of course much experience gained,
from tuning other systems, especially toy problems.

...

:- WHAT IS THE QUESTION AGAIN?
:  What is the most efficient and robust "skeletal" algorithm that allows
:  one to plug in any MATCH function, and it will find the global extremum.
:  I know many people will say it depends on the MATCH function, but I
:  would like to find out the central elements that make this TYPE of
:  optimization efficient and robust.

Yes, that is what most people will, and _should_ say.

...

:  The down side of this approach is that you now have more
:  variables to optimize, so your computing time is increased. Thus,
:  what is a good compromise???

This is extremely problem dependent.  Consider a simple object
moving in some n-dim quasi-paraboloid.  Newton's method in this
n-dim space will likely be more efficient than any sampling technique.
Now take an arbitrary projection into a m-dim space, m < n, and
the motion will look quite radical, and perhaps you will need a
sampling technique to find the global minimum.

Add noise to the data, and the problems really multiply.  I have
seen projects trying to optimize brain/head-shapes from MRI data,
with "standard" deterministic algorithms, fail to get the optimal
point 50% of the time.  Is this OK for your project?

More details of your "cost function" followed.  While quite relevant
to setting up your problem, if you are reasonably sure that this
is well formulated, then you likely will spend the next week or so
most efficiently by trying to work with an optimization program
per se.

:- HOW ABOUT QUATERNIONS? From my limited understanding, quaternions
:  allow you to specify rotation around arbitrary axis. However, I don't
:  know how it can represent an orientation ORI in a way different from
:  a 3x3 matrix. Any concrete suggestions? If there is a way, will it
:  overcome the problem of inter-dependence of rotations? Also how would
:  it handle translation, or do you just treat translation separately?

Quaternions are a 2-1 representation of the rotation group.  After
this constraint is put in, I don't see how this can help, but
perhaps this might lead to a more efficient way of actually writing
your cost function.

:- ANY OTHER WAYS OF ESTABLISHING INDEPENDENT VARIABLES FOR SPECIFYING
:  ORIENTATION?

Oh yes, oh yes.  See
	%A H. Goldstein
	%T Classical Mechanics, 2nd ed.
	%I Addison Wesley
	%C Reading, MA
	%D 1980

:II. OPTIMIZATION
:----------------
:
:- GENERAL DIFFICULTY of optimization for image coregistration is
:  that rotation causes non-linearity in the parameter space as
:  discussed above. Secondly, the match function usually has *MANY*
:  local maxima (for the MATCH function, or local minima for the DIFF
:  function), especially if the images are from pathological
:  subjects. Hence, there is a need to find robust algorithms for
:  avoiding these extrema.

If true, you likely need to use a global optimiztion code.  Most
likely, you need to have a code that will permit you to tune it to
your problem.  The downside of having a global optimization code
that can be tuned to be robust over a large class of problems, is
that the learning curve will be pretty steep on how to tune the
code and how to tune it for your problem.  How to best spend your
time is an administrative decision you must make based on the
quality of answer you need in how much time and resources you have
available.

:- PROBLEMS WITH GRADIENT METHODS. From personal experience (and I have
:  found many people to disagree with this), gradient methods
:  do not work well if the gradients cannot be calculated analytically.
:  In image coregistration, the derivative is (as far as I know)
:  approximated by finite different method by resampling the image at
:  slightly different orientations (denoted by x below):
:
:	dF/dx ~ (F(x+dx) - F(x))/dx
:
:  Many researchers have developed hacks to tune "dx" in the above
:  problem. If dx is fixed, I've found that gradient-descent methods often
:  get close to the global extremum, but not quite there. The quick
:  hacks for tuning dx usually work for a small group of well-defined
:  problems but not for anything else.

Generally, sampling codes do not require derivatives, although
some, like ASA, can optionally use them for some particular methods
to make the optimization more efficient.

:- ITERATIVE STEPS FOR ROTATION VERSUS TRANSLATION. When optimizing
:  orientation in 6 degrees of freedom, one inevitably asks the question,
:  should I rotate by X degrees or should I translate by Y mm when both
:  transformation will lead to a better match (i.e. higher value for the
:  MATCH function). THE QUESTION IS: what is the *functional* equivalence
:  between rotation and translation? i.e. If both the rotation and the
:  translation will improve the match by the SAME amount, which is the
:  better transformation to do? IS IT POSSIBLE TO DECIDE???

I cannot see how this cannot depend on the "landscape" of the cost function?

...

:- INITIAL ESTIMATES are often required to 1) reduce computing time
:  for optimization, 2) avoid local extrema. I have reported a fully
:  automated algorithm for determining brain orientation in 3D SPECT
:  without any initial estimates (J Nucl Med 37(5): 8P. 1996). However,
:  this algorithm is based on a priori knowledge (i.e. the rules programmed
:  by me) about the well conserved, overall geometry of the brain.
:  THE QUESTION IS: has anyone developed any generic algorithms for
:  resolving the initial pose (i.e. the approximate match between two
:  images.) Alternately, has anyone developed any algorithms fast enough
:  that it can be used to test ALL orientations?

This is a reasonable question that might be answered in the usenet
groups you have posted.  Has anyone solved/optimized at least parts
of your particular problem already?

:- STOCHASTIC OPTIMIZATION is very similar to genetic algorithm, but is not
:  as full-pledged as GA because it only alters the parameters in the
:  optimization algorithm, but not components of the algorithm itself.
:  Has anyone performed any comparisons between the two?

I think this is a fairy tale, at the least a very shallow comparison.
If you are already familiar with GA, then it makes sense to try
this approach first.  If you are already familiar with SA, then it
makes sense to try this approach first.  Comparisons between the
two are often with toy problems, as in
          %A L. Ingber
          %A B. Rosen
          %T Genetic algorithms and very fast simulated reannealing:
             A comparison
          %J Mathematical Computer Modelling
          %V 16
          %N 11
          %P 87-100
          %D 1992
          %O URL http://www.ingber.com/asa92_saga.ps.Z
or are embedded in specific contexts, as in
          %A D.G. Mayer
          %A J.A. Belward
          %A K. Burrage
          %A M.A Stuart
          %P 1995
          %T Optimization of a dairy farm model - Comparison of simulated
             annealing, simulated quenching and genetic algorithms
          %B Proceedings 1995 International Congress on Modelling and
          Simulation, 27-30 November 1995, University of Newcastle
          %P 33-38
          %D 1995

	  Contact David G. Mayer <davidm@ysun73.ari.dpi.qld.gov.au>
	  for further information.

:- ANY REAL DIFFERENCE AMONG 1) NEAREST NEIGHBOR, 2) TRILINEAR, 3) TRICUBIC,
:  AND 4) SINC INTERPOLATION? These 4 interpolation methods are listed in
:  the order of increasing accuracy. However, accuracy increases at the
:  expense of increased computation time. THE QUESTION IS: since one has
:  to sample many points for image coregistration, interpolation errors
:  would probably cancel out, so does it make sense to use the more accurate
:  interpolation methods?
:
:- DERIVATIVE CALCULATION. Suppose I want to find the derivative of the
:  MATCH function with respect to rotation around the X axis, the usual
:  way is to 1) rotate the image to the right a little and resample,
:	    2) rotate the image to the left  a little and resample,
:            3) calculate the MATCH function at each orientation,
:            4) and take the difference.
:  i.e. use the finite difference method. Is there a way of calculating
:  the derivative so that I don't have to do so many resampling
:  calculations?

Again, I think the effect on the cost function and its optimization
will be landscape-dependent.  However, it does seem that you should
be able to get answers from other sources on the inherent numerical
accuracies of these techniques.

:
:- ONLY 3 POINTS are needed (in most cases) to completely define an
:  object's orientation. Therefore, if we can confidently and accurately
:  define 3 corresponding points in each image (i.e. 3x2 = 6 points in
:  total), we can coregister the images already. In reality, we pick *a
:  lot* more points because medical images are noisy and we cannot
:  confidently localize corresponding positions on two different images.
:  The question is, how many points do we need for an *AUTOMATED* algorithm.
:  Surface-based registration algorithms show that you only need a subset
:  of all the points in a volume for coregistration, but numerous studies
:  have shown that surface-based methods are not as accurate as volume-
:  based methods, so obviously using the outermost surface does not seem
:  to be a good criterion for choosing points. This can be understood
:  intuitively via Pelizarri's hat-and-head analogy for surface fitting:
:  If the hat does not fit the head snugly (i.e. there is some gap between
:  2 surfaces to be matched), then one can turn or tilt the hat slightly
:  with noticing much difference in the match.
:  Thus, it seems more worthwhile to pick points where there are more
:  spatial variations. i.e. FEATURE POINTS. Definitions for feature points
:  abound. What is a good definition and a good algorithm for finding them?

Another slant on selecting surface versus volume optimization deals
with the parameterization of the problem.  If you actually are
trying to fit an m-dim surface, m < n, but are using an n dim volume
of data + adding constraints for the other n-m dims, then you should
try to explicity/numerically solve for the m-dim surface in your
cost function, so you wind up fitting only m parameters.  In the
context of a sampling technique, it is an extreme waste of sampling
to have to sample additional dimensions; most often it is a good
tradeoff to add som nonlinearlity by working with the m dim space.

:- TETRAHEDRAL LATTICE. Suppose that one abandons the use of feature points
:  as discussed in the previous paragraph and sample points on a regularly
:  structured lattice, one should realize that a rectangular grid (i.e.
:  the image matrix) is not the best way to sample a volume because some
:  unsampled positions are quite far from any of the lattice points. A
:  better "packing" of the sampled positions into a volume would be to
:  arrange them into a tetrahedral grid. Using a tetrahedral lattice is
:  a nuisance for programming, but it does seem cover a much better sampling
:  volume. Any comments on this?

It seems your main concern is CPU time for optimization?  If so,
the labor of programming is a relatively small concern.  If not,
then there are administrative decisions to be made similar to those
mentioned above.

:- MULTI-RESOLUTION is useful for reducing computation time. I haven't found
:  it too useful for avoiding local extrema. I've never seen it done so
:  that, for example, if you rotate an image 180 degs, multi-resolution
:  will get you out of that local extrema. Does anyone have any other
:  opinions?

This depends on how much you know about your system.  For example,
if you have some good idea of the relative values of hills and
valleys, then you could do a first optimization pass using a coarse
resolution ("coarse" specific to your problem), and then follow up
with finer global or possibly even quasi-local optimization.  I
have done this in several systems, including EEG analyses,
          %A L. Ingber
          %T Statistical mechanics of neocortical interactions:
             Canonical momenta indicators of electroencephalography
          %J Physical Review E
          %P (to be published)
          %D 1997
          %O URL http://www.ingber.com/smni97_cmi.ps.Z
financial markets,
          %A L. Ingber
          %T Canonical momenta indicators of financial markets and
             neocortical EEG
          %B International Conference on Neural Information Processing
             (ICONIP'96)
          %I Springer
          %C New York
          %P 777-784
          %D 1996
          %O Invited paper to the 1996 International Conference on Neural
             Information Processing (ICONIP'96), Hong Kong, 24-27 September
             1996. URL http://www.ingber.com/markets96_momenta.ps.Z
and combat analysis
          %A M. Bowman
          %A L. Ingber
          %T Canonical momenta of nonlinear combat
          %B Proceedings of the 1997 Simulation Multi-Conference, 6-10
             April 1997, Atlanta, GA
          %I Society for Computer Simulation
          %C San Diego, CA
          %P (to be published)
          %D 1997
          %O URL http://www.ingber.com/combat97_cmi.ps.Z


:---------------------------------------------------------------------------
:
:Thank you very much for reading this ***long*** posting. Looking forward
:to your responses.

And we look forward to your ***long*** summary of replies :>).

Lester

-- 
/*             RESEARCH                            ingber@ingber.com *
 *       INGBER                                 ftp://ftp.ingber.com *
 * LESTER                                     http://www.ingber.com/ *
 * Prof. Lester Ingber __ PO Box 857 __ McLean, VA 22101-0857 __ USA */