implement optimization using newton method

implement optimization using newton method

6 files changed, 132 insertions, 11 deletions

diff --git a/src/GipfelWidget.H b/src/GipfelWidget.H
index cd803d4..2d05f2c 100644
--- a/src/GipfelWidget.H
+++ b/src/GipfelWidget.H
@@ -1,5 +1,5 @@
 // 
-// "$Id: GipfelWidget.H,v 1.8 2005/04/19 16:17:56 hofmann Exp $"
+// "$Id: GipfelWidget.H,v 1.9 2005/04/24 09:51:17 hofmann Exp $"
 //
 // X11 header file for the Fast Light Tool Kit (FLTK).
 //
@@ -60,6 +60,8 @@ class GipfelWidget : public Fl_Widget {
   
   int comp_params();
 
+  int newton();
+
   void draw();
 };
 #endif
diff --git a/src/GipfelWidget.cxx b/src/GipfelWidget.cxx
index b43ea86..365def1 100644
--- a/src/GipfelWidget.cxx
+++ b/src/GipfelWidget.cxx
@@ -1,5 +1,5 @@
 // 
-// "$Id: GipfelWidget.cxx,v 1.10 2005/04/19 16:17:56 hofmann Exp $"
+// "$Id: GipfelWidget.cxx,v 1.11 2005/04/24 09:51:17 hofmann Exp $"
 //
 // PSEditWidget routines.
 //
@@ -178,6 +178,13 @@ GipfelWidget::comp_params() {
   redraw();
 }
 
+
+int
+GipfelWidget::newton() {
+  pan->newton();
+  redraw();
+}
+
 int
 GipfelWidget::handle(int event) {
   int mark_x, mark_y;
diff --git a/src/Makefile b/src/Makefile
index e42ae7b..79484f6 100644
--- a/src/Makefile
+++ b/src/Makefile
@@ -1,7 +1,7 @@
 CC=gcc
 CPP=c++
 CFLAGS=-g -I/usr/X11R6/include
-CXXFLAGS=-g -I/usr/X11R6/include -fno-stack-protector
+CXXFLAGS=-g -I/usr/X11R6/include -I/usr/local/include -fno-stack-protector
 LDFLAGS=-g -fno-stack-protector
 
 OBJECTS=flmountains.o GipfelWidget.o Mountain.o Panorama.o
@@ -13,7 +13,7 @@ OBJECTS=flmountains.o GipfelWidget.o Mountain.o Panorama.o
 	$(CPP) -c $(CXXFLAGS) $*.cxx
 
 gipfel: $(OBJECTS)
-	$(CC) -fno-stack-protector -o gipfel $(OBJECTS) -L/usr/X11R6/lib -lm -lfltk -lfltk_images
+	$(CC) -fno-stack-protector -o gipfel $(OBJECTS) -L/usr/X11R6/lib -L/usr/local/lib -lm -lccm -lfltk -lfltk_images
 
 flmountains: $(flmountains.o)
 	$(CPP) -o flmountains flmountains.o -L/usr/X11R6/lib -lm -lfltk -lfltk_images
diff --git a/src/Panorama.H b/src/Panorama.H
index 12586a8..24392ec 100644
--- a/src/Panorama.H
+++ b/src/Panorama.H
@@ -1,5 +1,5 @@
 // 
-// "$Id: Panorama.H,v 1.6 2005/04/19 16:17:56 hofmann Exp $"
+// "$Id: Panorama.H,v 1.7 2005/04/24 09:51:17 hofmann Exp $"
 //
 // X11 header file for the Fast Light Tool Kit (FLTK).
 //
@@ -38,7 +38,11 @@ class Panorama {
   double scale;
   double a_nick;
   double a_tilt;
-  
+  // these must be set before using newton()
+  double tan_dir_view, tan_nick_view, n_scale;
+  double tan_dir_m1, tan_nick_m1, tan_dir_m2, tan_nick_m2;
+  double d_m1_2, d_m2_2, d_m1_m2_2;
+
   int get_pos(const char *name, double *phi, double *lam, double *height);
 
   void update_visible_mountains();
@@ -57,6 +61,10 @@ class Panorama {
 
   double comp_scale(double alph_a, double alph_b, double d1, double d2);
 
+  int get_matrix(double m[]);
+
+
+
  public:
   Panorama();
 
@@ -81,5 +89,7 @@ class Panorama {
   void set_scale(double s);
 
   int comp_params();
+
+  int newton();
 };
 #endif
diff --git a/src/Panorama.cxx b/src/Panorama.cxx
index f812e61..f68132b 100644
--- a/src/Panorama.cxx
+++ b/src/Panorama.cxx
@@ -1,5 +1,5 @@
 // 
-// "$Id: Panorama.cxx,v 1.8 2005/04/19 16:17:56 hofmann Exp $"
+// "$Id: Panorama.cxx,v 1.9 2005/04/24 09:51:17 hofmann Exp $"
 //
 // PSEditWidget routines.
 //
@@ -25,7 +25,9 @@
 #include <stdlib.h>
 #include <string.h>
 #include <math.h>
-
+extern "C" {
+#include <ccmath.h>
+}
 #include "Panorama.H"
 
 Panorama::Panorama() {
@@ -114,7 +116,9 @@ Panorama::get_visible_mountains() {
 int
 Panorama::set_mountain(Mountain *m, int x, int y) {
   if (m1 && m2 && m != m1 && m != m2) {
-    return 1;
+    fprintf(stderr, "Resetting mountains\n");
+    m1 = NULL;
+    m2 = NULL;
   }
 
   m->x = x;
@@ -147,9 +151,10 @@ Panorama::comp_params() {
   //  a_nick = comp_center_angle(m1->a_nick, m2->a_nick, -m1->y, -m2->y) - pi;  
   fprintf(stderr, "center = %f, scale = %f, nick=%f\n", a_center /deg2rad, scale, a_nick/deg2rad);
   update_visible_mountains();
-
+#if 0
   m1 = NULL;
   m2 = NULL;
+#endif
   return 0;
 }
 
@@ -345,3 +350,92 @@ Panorama::nick(double dist, double height) {
   
   return beta - pi / 2.0;
 }
+
+
+
+
+
+//
+//
+//
+
+int
+Panorama::get_matrix(double m[]) {
+  
+  m[0] = pow(n_scale,2.0)*(1.0/pow((tan_nick_m1*tan_nick_view + 1.0),2.0)*(2.0*tan_nick_m1 - 2.0 * tan_nick_view) + 2.0*tan_nick_m1*pow((tan_nick_m1 - tan_nick_view), 2.0)/pow((tan_nick_m1*tan_nick_view + 1.0), 3.0));
+
+  m[3] = pow(n_scale, 2.0) *(1.0/pow((tan_dir_m1*tan_dir_view + 1.0), 2.0) * (2.0*tan_dir_m1 - 2.0*tan_dir_view) + 2.0*tan_dir_m1*pow((tan_dir_m1 - tan_dir_view),2.0) / pow((tan_dir_m1*tan_dir_view + 1.0), 3.0));
+
+  m[6] = -2.0*n_scale*(pow((tan_dir_m1 - tan_dir_view), 2.0)/pow((tan_dir_m1*tan_dir_view + 1.0), 2.0) + pow((tan_nick_m1 - tan_nick_view), 2.0)/pow((tan_nick_m1*tan_nick_view + 1.0), 2.0));
+
+  m[1] = pow(n_scale, 2.0)*(1.0/pow((tan_nick_m2*tan_nick_view + 1.0), 2.0)*(2.0*tan_nick_m2 - 2.0*tan_nick_view) + 2.0*tan_nick_m2*pow((tan_nick_m2 - tan_nick_view), 2.0)/pow((tan_nick_m2*tan_nick_view + 1.0), 3.0));
+
+  m[4] = pow(n_scale, 2.0)*(1.0/pow((tan_dir_m2*tan_dir_view + 1.0), 2.0)*(2.0*tan_dir_m2 - 2.0*tan_dir_view) + 2.0*tan_dir_m2*pow((tan_dir_m2 - tan_dir_view), 2.0)/pow((tan_dir_m2*tan_dir_view + 1.0), 3.0));
+
+  m[7] = -2.0*n_scale*(pow((tan_dir_m2 - tan_dir_view), 2.0)/pow((tan_dir_m2*tan_dir_view + 1.0), 2.0) + pow((tan_nick_m2 - tan_nick_view), 2.0)/pow((tan_nick_m2*tan_nick_view + 1.0), 2.0));
+
+  m[2] = 2.0*(n_scale*(tan_nick_m1 - tan_nick_view)/(tan_nick_m1*tan_nick_view + 1.0) - n_scale*(tan_nick_m2 - tan_nick_view)/(tan_nick_m2*tan_nick_view + 1.0))*(n_scale/(tan_nick_m1*tan_nick_view + 1.0) - n_scale/(tan_nick_m2*tan_nick_view + 1.0) + tan_nick_m1*n_scale*(tan_nick_m1 - tan_nick_view)/pow((tan_nick_m1*tan_nick_view + 1.0), 2.0) - tan_nick_m2*n_scale*(tan_nick_m2 - tan_nick_view)/pow((tan_nick_m2*tan_nick_view + 1.0),2.0));
+
+  m[5] = 2.0*(n_scale*(tan_dir_m1 - tan_dir_view)/(tan_dir_m1*tan_dir_view + 1.0) - n_scale*(tan_dir_m2 - tan_dir_view)/(tan_dir_m2*tan_dir_view + 1.0))*(n_scale/(tan_dir_m1*tan_dir_view + 1.0) - n_scale/(tan_dir_m2*tan_dir_view + 1.0) + tan_dir_m1*n_scale*(tan_dir_m1 - tan_dir_view)/pow((tan_dir_m1*tan_dir_view + 1.0), 2.0) - tan_dir_m2*n_scale*(tan_dir_m2 - tan_dir_view)/pow((tan_dir_m2*tan_dir_view + 1.0), 2.0));
+
+  m[8] = - 2.0*(n_scale*(tan_dir_m1 - tan_dir_view)/(tan_dir_m1*tan_dir_view + 1.0) - n_scale*(tan_dir_m2 - tan_dir_view)/(tan_dir_m2*tan_dir_view + 1.0))*((tan_dir_m1 - tan_dir_view)/(tan_dir_m1*tan_dir_view + 1.0) - (tan_dir_m2 - tan_dir_view)/(tan_dir_m2*tan_dir_view + 1.0)) - 2.0*(n_scale*(tan_nick_m1 - tan_nick_view)/(tan_nick_m1*tan_nick_view + 1.0) - n_scale*(tan_nick_m2 - tan_nick_view)/(tan_nick_m2*tan_nick_view + 1.0))*((tan_nick_m1 - tan_nick_view)/(tan_nick_m1*tan_nick_view + 1.0) - (tan_nick_m2 - tan_nick_view)/(tan_nick_m2*tan_nick_view + 1.0));
+
+  return 0;
+}
+
+
+int 
+Panorama::newton() {
+  double a[9];
+  double b[3];
+  double a_x0[3], f_x0 [3], x0[3];
+  int ret;
+
+  tan_nick_view = tan(a_nick);
+  tan_dir_view = tan(a_center);
+  n_scale = scale;
+  tan_dir_m1 = tan(m1->alph);
+  tan_nick_m1 = tan(m1->a_nick);
+  tan_dir_m2 = tan(m2->alph);
+  tan_nick_m2 = tan(m2->a_nick);
+
+  fprintf(stderr, "m1: %d, %d; m2: %d, %d\n", m1->x, m1->y, m2->x, m2->y);
+  d_m1_2 = pow(m1->x, 2.0) + pow(m1->y, 2.0);
+  d_m2_2 = pow(m2->x, 2.0) + pow(m2->y, 2.0);
+  d_m1_m2_2 = pow(m1->x - m2->x, 2.0) + pow(m1->y - m2->y, 2.0);
+
+  fprintf(stderr, "d_m1_2 %f, d_m2_2 %f, d_m1_m2_2 %f\n", 
+	  d_m1_2, d_m2_2, d_m1_m2_2);
+
+  get_matrix(a);
+
+  f_x0[0] = d_m1_2 - (pow((tan_nick_view-tan_nick_m1),2.0)/pow((tan_nick_m1*tan_nick_view+1), 2.0)+pow((tan_dir_view-tan_dir_m1),2.0)/pow((tan_dir_m1*tan_dir_view+1),2.0))*pow(n_scale, 2.0);
+
+  f_x0[1] = d_m2_2 - (pow((tan_nick_view-tan_nick_m2),2.0)/pow((tan_nick_m2*tan_nick_view+1),2.0)+pow((tan_dir_view-tan_dir_m2),2.0)/pow((tan_dir_m2*tan_dir_view+1),2.0))*pow(n_scale, 2.0);
+
+  f_x0[2] = d_m1_m2_2 - (pow((- (((tan_dir_view - tan_dir_m1) * n_scale) / (tan_dir_m1 * tan_dir_view + 1.0)) + (((tan_dir_view - tan_dir_m2) * n_scale) / (tan_dir_m2 * tan_dir_view + 1))), 2.0) + pow((- (((tan_nick_view - tan_nick_m1) * n_scale) / (tan_nick_m1 * tan_nick_view + 1)) + ((tan_nick_view - tan_nick_m2) * n_scale) / (tan_nick_m2 * tan_nick_view + 1)), 2.0));
+
+  fprintf(stderr, "f_x0[0] %f, f_x0[1] %f, f_x0[2] %f\n", 
+	  f_x0[0], f_x0[1], f_x0[2]);
+
+  x0[0] = tan_nick_view;
+  x0[1] = tan_dir_view;
+  x0[2] = n_scale;
+
+  rmmult(a_x0, a, x0, 3, 3, 1);
+    
+  b[0] = a_x0[0] - f_x0[0];
+  b[1] = a_x0[1] - f_x0[1];
+  b[2] = a_x0[2] - f_x0[2];
+
+  
+  ret = solv(a, b, 3);
+  fprintf(stderr, "solv returned %d\n", ret);
+
+  a_nick = atan(b[0]);
+  a_center = atan(b[1]);
+  scale = b[2];
+  fprintf(stderr, "center = %f, scale = %f, nick=%f\n", a_center /deg2rad, scale, a_nick/deg2rad);
+  update_visible_mountains();
+  return 0;
+}
diff --git a/src/gipfel.cxx b/src/gipfel.cxx
index 934ed75..8335abb 100644
--- a/src/gipfel.cxx
+++ b/src/gipfel.cxx
@@ -1,5 +1,5 @@
 // 
-// "$Id: gipfel.cxx,v 1.11 2005/04/19 16:17:56 hofmann Exp $"
+// "$Id: gipfel.cxx,v 1.12 2005/04/24 09:51:17 hofmann Exp $"
 //
 // flpsed program.
 //
@@ -90,6 +90,12 @@ void comp_cb(Fl_Widget *, void *) {
   }
 }
 
+void newton_cb(Fl_Widget *, void *) {
+  if (gipf) {
+    gipf->newton();
+  }
+}
+
 void about_cb() {
   fl_message("flpsed -- a pseudo PostScript editor\n"
 	     "(c) Johannes Hofmann 2004, 2005\n\n"
@@ -207,6 +213,8 @@ int main(int argc, char** argv) {
   t->align(FL_ALIGN_LEFT);
   Fl_Button *b = new Fl_Button(200, 60, 20, 15, "comp");
   b->callback(comp_cb);
+  Fl_Button *b1 = new Fl_Button(250, 60, 20, 15, "opt");
+  b1->callback(newton_cb);
 
   scroll = new Fl_Scroll(0, 75, win->w(), win->h()-75);