--- /dev/null
+/***********************************************************************
+ * Software License Agreement (BSD License)
+ *
+ * Copyright 2008-2009 Marius Muja (mariusm@cs.ubc.ca). All rights reserved.
+ * Copyright 2008-2009 David G. Lowe (lowe@cs.ubc.ca). All rights reserved.
+ *
+ * THE BSD LICENSE
+ *
+ * Redistribution and use in source and binary forms, with or without
+ * modification, are permitted provided that the following conditions
+ * are met:
+ *
+ * 1. Redistributions of source code must retain the above copyright
+ * notice, this list of conditions and the following disclaimer.
+ * 2. Redistributions in binary form must reproduce the above copyright
+ * notice, this list of conditions and the following disclaimer in the
+ * documentation and/or other materials provided with the distribution.
+ *
+ * THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR
+ * IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
+ * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
+ * IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY DIRECT, INDIRECT,
+ * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT
+ * NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+ * DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+ * THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+ * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF
+ * THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+ *************************************************************************/
+
+#ifndef SIMPLEX_DOWNHILL_H
+#define SIMPLEX_DOWNHILL_H
+
+namespace flann
+{
+
+/**
+ Adds val to array vals (and point to array points) and keeping the arrays sorted by vals.
+*/
+template <typename T>
+void addValue(int pos, float val, float* vals, T* point, T* points, int n)
+{
+ vals[pos] = val;
+ for (int i=0;i<n;++i) {
+ points[pos*n+i] = point[i];
+ }
+
+ // bubble down
+ int j=pos;
+ while (j>0 && vals[j]<vals[j-1]) {
+ swap(vals[j],vals[j-1]);
+ for (int i=0;i<n;++i) {
+ swap(points[j*n+i],points[(j-1)*n+i]);
+ }
+ --j;
+ }
+}
+
+
+/**
+ Simplex downhill optimization function.
+ Preconditions: points is a 2D mattrix of size (n+1) x n
+ func is the cost function taking n an array of n params and returning float
+ vals is the cost function in the n+1 simplex points, if NULL it will be computed
+
+ Postcondition: returns optimum value and points[0..n] are the optimum parameters
+*/
+template <typename T, typename F>
+float optimizeSimplexDownhill(T* points, int n, F func, float* vals = NULL )
+{
+ const int MAX_ITERATIONS = 10;
+
+ assert(n>0);
+
+ T* p_o = new T[n];
+ T* p_r = new T[n];
+ T* p_e = new T[n];
+
+ int alpha = 1;
+
+ int iterations = 0;
+
+ bool ownVals = false;
+ if (vals == NULL) {
+ ownVals = true;
+ vals = new float[n+1];
+ for (int i=0;i<n+1;++i) {
+ float val = func(points+i*n);
+ addValue(i, val, vals, points+i*n, points, n);
+ }
+ }
+ int nn = n*n;
+
+ while (true) {
+
+ if (iterations++ > MAX_ITERATIONS) break;
+
+ // compute average of simplex points (except the highest point)
+ for (int j=0;j<n;++j) {
+ p_o[j] = 0;
+ for (int i=0;i<n;++i) {
+ p_o[i] += points[j*n+i];
+ }
+ }
+ for (int i=0;i<n;++i) {
+ p_o[i] /= n;
+ }
+
+ bool converged = true;
+ for (int i=0;i<n;++i) {
+ if (p_o[i] != points[nn+i]) {
+ converged = false;
+ }
+ }
+ if (converged) break;
+
+ // trying a reflection
+ for (int i=0;i<n;++i) {
+ p_r[i] = p_o[i] + alpha*(p_o[i]-points[nn+i]);
+ }
+ float val_r = func(p_r);
+
+ if (val_r>=vals[0] && val_r<vals[n]) {
+ // reflection between second highest and lowest
+ // add it to the simplex
+ logger.info("Choosing reflection\n");
+ addValue(n, val_r,vals, p_r, points, n);
+ continue;
+ }
+
+ if (val_r<vals[0]) {
+ // value is smaller than smalest in simplex
+
+ // expand some more to see if it drops further
+ for (int i=0;i<n;++i) {
+ p_e[i] = 2*p_r[i]-p_o[i];
+ }
+ float val_e = func(p_e);
+
+ if (val_e<val_r) {
+ logger.info("Choosing reflection and expansion\n");
+ addValue(n, val_e,vals,p_e,points,n);
+ }
+ else {
+ logger.info("Choosing reflection\n");
+ addValue(n, val_r,vals,p_r,points,n);
+ }
+ continue;
+ }
+ if (val_r>=vals[n]) {
+ for (int i=0;i<n;++i) {
+ p_e[i] = (p_o[i]+points[nn+i])/2;
+ }
+ float val_e = func(p_e);
+
+ if (val_e<vals[n]) {
+ logger.info("Choosing contraction\n");
+ addValue(n,val_e,vals,p_e,points,n);
+ continue;
+ }
+ }
+ {
+ logger.info("Full contraction\n");
+ for (int j=1;j<=n;++j) {
+ for (int i=0;i<n;++i) {
+ points[j*n+i] = (points[j*n+i]+points[i])/2;
+ }
+ float val = func(points+j*n);
+ addValue(j,val,vals,points+j*n,points,n);
+ }
+ }
+ }
+
+ float bestVal = vals[0];
+
+ delete[] p_r;
+ delete[] p_o;
+ delete[] p_e;
+ if (ownVals) delete[] vals;
+
+ return bestVal;
+}
+
+}
+
+#endif //SIMPLEX_DOWNHILL_H
projects / opencv / blobdiff