summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--configure.ac9
-rw-r--r--src/GipfelWidget.H4
-rw-r--r--src/GipfelWidget.cxx122
-rw-r--r--src/Hill.H12
-rw-r--r--src/Hill.cxx27
-rw-r--r--src/ImageMetaData.H6
-rw-r--r--src/ImageMetaData.cxx34
-rw-r--r--src/Makefile.am2
-rw-r--r--src/Makefile.lsq_funcs4
-rw-r--r--src/Panorama.H5
-rw-r--r--src/Panorama.cxx23
-rw-r--r--src/Projection.H2
-rw-r--r--src/Projection.cxx2
-rw-r--r--src/ProjectionSphaeric.H2
-rw-r--r--src/ProjectionSphaeric.cxx12
-rw-r--r--src/ProjectionTangential.H2
-rw-r--r--src/ProjectionTangential.cxx12
-rw-r--r--src/ProjectionTangentialLSQ.H31
-rw-r--r--src/ProjectionTangentialLSQ.cxx355
-rw-r--r--src/ViewParams.H2
-rw-r--r--src/lsq_funcs.mac40
21 files changed, 615 insertions, 93 deletions
diff --git a/configure.ac b/configure.ac
index 27e4908..347ff2b 100644
--- a/configure.ac
+++ b/configure.ac
@@ -47,6 +47,15 @@ fi
CXXFLAGS="`$FLTKCONFIG --use-images --cflags` $CXXFLAGS"
LIBS="`$FLTKCONFIG --use-images --ldflags` $LIBS"
+# Check for gsl
+AC_PATH_PROG(GSLCONFIG,gsl-config)
+if test "x$GSLCONFIG" = x; then
+ echo "gsl-config not found"
+ exit 1
+fi
+
+CXXFLAGS="`$GSLCONFIG --cflags` $CXXFLAGS"
+LIBS="`$GSLCONFIG --libs` $LIBS"
# Check for ccmath
AC_CHECK_HEADERS([ccmath.h], [], [echo "Error: ccmath.h not found."; exit 1;])
diff --git a/src/GipfelWidget.H b/src/GipfelWidget.H
index 80761a5..cba4513 100644
--- a/src/GipfelWidget.H
+++ b/src/GipfelWidget.H
@@ -17,7 +17,7 @@ class GipfelWidget : public Fl_Widget {
Hill *cur_mountain;
Hills *marker;
Hills *track_points;
- Hill *m1, *m2;
+ Hills *known_hills;
Panorama *pan;
Fl_Menu_Button *mb;
char *img_file;
@@ -28,6 +28,8 @@ class GipfelWidget : public Fl_Widget {
int set_cur_mountain(int m_x, int m_y);
+ int toggle_known_mountain(int m_x, int m_y);
+
int set_mountain(int m_x, int m_y);
void set_labels(Hills *v);
diff --git a/src/GipfelWidget.cxx b/src/GipfelWidget.cxx
index 00c5aad..e1dad77 100644
--- a/src/GipfelWidget.cxx
+++ b/src/GipfelWidget.cxx
@@ -46,8 +46,7 @@ GipfelWidget::GipfelWidget(int X,int Y,int W, int H): Fl_Widget(X, Y, W, H) {
cur_mountain = NULL;
mb = NULL;
marker = new Hills();
- m1 = NULL;
- m2 = NULL;
+ known_hills = new Hills();
img_file = NULL;
track_width = 200.0;
show_hidden = 0;
@@ -102,7 +101,8 @@ GipfelWidget::load_image(char *file) {
set_tilt_angle(md->get_tilt());
set_projection((Projection::Projection_t) md->get_projection_type());
set_focal_length_35mm(md->get_focal_length_35mm());
-
+ md->get_distortion_params(&pan->parms.k0, &pan->parms.k1);
+ fprintf(stderr, "%f %f\n", pan->parms.k0, pan->parms.k1);
delete md;
return 0;
@@ -133,6 +133,7 @@ GipfelWidget::save_image(char *file) {
md->set_tilt(get_tilt_angle());
md->set_focal_length_35mm(get_focal_length_35mm());
md->set_projection_type((int) get_projection());
+ md->set_distortion_params(pan->parms.k0, pan->parms.k1);
ret = md->save_image(img_file, file);
delete md;
@@ -211,8 +212,6 @@ void
GipfelWidget::draw() {
Hills *mnts;
Hill *m;
- int center_x = w() / 2;
- int center_y = h() / 2;
int i;
if (img == NULL) {
@@ -227,6 +226,8 @@ GipfelWidget::draw() {
mnts = pan->get_visible_mountains();
for (i=0; i<mnts->get_num(); i++) {
m = mnts->get(i);
+ int m_x = w() / 2 + x() + (int) rint(m->x);
+ int m_y = h() / 2 + y() + (int) rint(m->y);
if (m->flags & (Hill::DUPLICATE|Hill::TRACK_POINT)) {
continue;
@@ -236,34 +237,19 @@ GipfelWidget::draw() {
continue;
}
- if (m == m1) {
+ if (known_hills->contains(m)) {
fl_color(FL_RED);
- draw_flag(center_x + m->x + x(), center_y + m->y + y(), "1");
- } else if (m == m2) {
- fl_color(FL_RED);
- draw_flag(center_x + m->x + x(), center_y + m->y + y(), "2");
+ draw_flag(m_x, m_y, "1");
} else if (m->flags & Hill::HIDDEN) {
fl_color(FL_BLUE);
} else {
fl_color(FL_BLACK);
}
- fl_xyline(center_x + m->x + x() - CROSS_SIZE, center_y + m->y + y(), center_x + m->x + x() + CROSS_SIZE);
- fl_yxline(center_x + m->x + x(), center_y + m->y + m->label_y + y() - CROSS_SIZE, center_y + m->y + y() + CROSS_SIZE);
-
- fl_draw(m->name,
- center_x + m->x + x(),
- center_y + m->y + m->label_y + y());
- }
-
- /* markers */
- for (i=0; i<marker->get_num(); i++) {
- m = marker->get(i);
+ fl_xyline(m_x- CROSS_SIZE, m_y, m_x + CROSS_SIZE);
+ fl_yxline(m_x, m_y + m->label_y - CROSS_SIZE, m_y + CROSS_SIZE);
- fl_color(FL_GREEN);
- fl_xyline(center_x + m->x + x() - CROSS_SIZE * 2, center_y + m->y + y(), center_x + m->x + x() + CROSS_SIZE * 2);
- fl_yxline(center_x + m->x + x(), center_y + m->y + y() - CROSS_SIZE * 2, center_y + m->y + y() + CROSS_SIZE * 2);
- draw_flag(center_x + m->x + x(), center_y + m->y + y(), NULL);
+ fl_draw(m->name, m_x , m_y + m->label_y);
}
/* track */
@@ -271,6 +257,9 @@ GipfelWidget::draw() {
int last_x, last_y, last_initialized = 0;
for (i=1; i<track_points->get_num(); i++) {
+ m = mnts->get(i);
+ int m_x = w() / 2 + x() + (int) rint(m->x);
+ int m_y = h() / 2 + y() + (int) rint(m->y);
if (!(track_points->get(i)->flags & Hill::VISIBLE)) {
continue;
}
@@ -285,14 +274,13 @@ GipfelWidget::draw() {
get_rel_track_width(track_points->get(i)));
if (last_initialized) {
fl_begin_line();
- fl_vertex(center_x + x() + last_x, center_y + y() + last_y);
- fl_vertex(center_x + x() + track_points->get(i)->x,
- center_y + y() + track_points->get(i)->y);
+ fl_vertex(last_x, last_y);
+ fl_vertex(m_x, m_y);
fl_end_line();
}
- last_x = track_points->get(i)->x;
- last_y = track_points->get(i)->y;
+ last_x = m_x;
+ last_y = m_y;
last_initialized++;
}
fl_line_style(0);
@@ -303,7 +291,7 @@ GipfelWidget::draw() {
static int
-overlap(int m1, int n1, int m2, int n2) {
+overlap(double m1, double n1, double m2, double n2) {
return m1 <= n2 && n1 >= m2;
}
@@ -354,6 +342,28 @@ GipfelWidget::set_labels(Hills *v) {
int
GipfelWidget::set_cur_mountain(int m_x, int m_y) {
+ Hill *m;
+ int center_x = w() / 2;
+ int center_y = h() / 2;
+ int i;
+
+ for (i=0; i<known_hills->get_num(); i++) {
+ m = known_hills->get(i);
+
+ if (m_x - center_x >= m->x - 2 && m_x - center_x < m->x + 2 &&
+ m_y - center_y >= m->y - 2 && m_y - center_y < m->y + 2) {
+
+ cur_mountain = m;
+ return 0;
+ }
+ }
+
+ return 1;
+}
+
+
+int
+GipfelWidget::toggle_known_mountain(int m_x, int m_y) {
Hills *mnts = pan->get_visible_mountains();
Hill *m;
int center_x = w() / 2;
@@ -368,19 +378,12 @@ GipfelWidget::set_cur_mountain(int m_x, int m_y) {
if (m_x - center_x >= m->x - 2 && m_x - center_x < m->x + 2 &&
m_y - center_y >= m->y - 2 && m_y - center_y < m->y + 2) {
- cur_mountain = m;
- if (m1 != NULL && m2 != NULL) {
- fprintf(stderr, "Resetting m1 and m2\n");
- m1 = NULL;
- m2 = NULL;
- }
- if (m1 == NULL) {
- m1 = cur_mountain;
- fprintf(stderr, "m1 = %s\n", m1->name);
- } else if (m2 == NULL) {
- m2 = cur_mountain;
- fprintf(stderr, "m2 = %s\n", m2->name);
+
+ if (known_hills->contains(m)) {
+ known_hills->remove(m);
+ } else {
+ known_hills->add(m);
}
redraw();
@@ -524,11 +527,11 @@ GipfelWidget::center() {
Hill *m = choose_hill(pan->get_close_mountains(), "Center Peak");
if (m) {
set_center_angle(m->alph / deg2rad);
- if (!m1 || (m1 && m2)) {
- m1 = m;
- } else {
- m2 = m;
+
+ if (!known_hills->contains(m)) {
+ known_hills->add(m);
}
+
}
}
@@ -589,12 +592,12 @@ GipfelWidget::get_mountains() {
int
GipfelWidget::comp_params() {
- if (m1 == NULL || m2 == NULL) {
+ if (known_hills->get_num() < 2) {
fprintf(stderr, "Position m1 and m2 first.\n");
return 1;
}
fl_cursor(FL_CURSOR_WAIT);
- pan->comp_params(m1, m2);
+ pan->comp_params(known_hills);
set_labels(pan->get_visible_mountains());
redraw();
fl_cursor(FL_CURSOR_DEFAULT);
@@ -602,12 +605,12 @@ GipfelWidget::comp_params() {
int
GipfelWidget::guess() {
- if (m1 == NULL) {
+ if (known_hills->get_num() < 1) {
fprintf(stderr, "Position m1 first.\n");
return 1;
}
fl_cursor(FL_CURSOR_WAIT);
- pan->guess(marker, m1);
+ pan->guess(marker, known_hills->get(0));
set_labels(pan->get_visible_mountains());
redraw();
fl_cursor(FL_CURSOR_DEFAULT);
@@ -639,20 +642,25 @@ GipfelWidget::handle(int event) {
switch(event) {
case FL_PUSH:
+ mark_x = Fl::event_x()-x();
+ mark_y = Fl::event_y()-y();
if (Fl::event_button() == 1) {
-
- mark_x = Fl::event_x()-x();
- mark_y = Fl::event_y()-y();
set_cur_mountain(mark_x, mark_y);
-
- Fl::focus(this);
- return 1;
+ } else if (Fl::event_button() == 2) {
+ toggle_known_mountain(mark_x, mark_y);
}
+
+ Fl::focus(this);
+ return 1;
break;
case FL_DRAG:
set_mountain(Fl::event_x()-x(), Fl::event_y()-y());
return 1;
break;
+ case FL_RELEASE:
+ cur_mountain = NULL;
+ return 1;
+ break;
case FL_FOCUS:
return 1;
break;
diff --git a/src/Hill.H b/src/Hill.H
index 2e95def..b181622 100644
--- a/src/Hill.H
+++ b/src/Hill.H
@@ -27,7 +27,7 @@ class Hill {
double a_view;
double a_nick;
double dist;
- int x, y;
+ double x, y;
int label_x, label_y;
char *name;
int flags;
@@ -36,7 +36,7 @@ class Hill {
Hill(const Hill& h);
- Hill(int x_tmp, int y_tmp);
+ Hill(double x_tmp, double y_tmp);
~Hill();
};
@@ -60,6 +60,8 @@ class Hills {
void add(Hill *m);
+ void remove(const Hill *m);
+
void add(Hills *h);
void sort_phi();
@@ -71,9 +73,11 @@ class Hills {
void clear();
void clobber();
+
+ int contains(const Hill *m) const;
- int get_num();
+ int get_num() const;
- Hill *get(int n);
+ Hill *get(int n) const;
};
#endif
diff --git a/src/Hill.cxx b/src/Hill.cxx
index d5ce1f6..0504eeb 100644
--- a/src/Hill.cxx
+++ b/src/Hill.cxx
@@ -40,7 +40,7 @@ Hill::Hill(const Hill& h) {
flags = h.flags;
}
-Hill::Hill(int x_tmp, int y_tmp) {
+Hill::Hill(double x_tmp, double y_tmp) {
name = NULL;
phi = 0.0;
lam = 0.0;
@@ -269,6 +269,27 @@ Hills::clear() {
num = 0;
}
+int
+Hills::contains(const Hill *m) const {
+ for(int i=0; i<get_num();i++) {
+ if (get(i) == m) {
+ return 1;
+ }
+ }
+
+ return 0;
+}
+
+void
+Hills::remove(const Hill *h) {
+ for(int i=0; i<get_num();i++) {
+ if (get(i) == h) {
+ memmove(&m[i], &m[i+1], (get_num() - i - 1) * sizeof(Hill*));
+ num--;
+ }
+ }
+}
+
void
Hills::clobber() {
int i;
@@ -283,12 +304,12 @@ Hills::clobber() {
}
int
-Hills::get_num() {
+Hills::get_num() const {
return num;
}
Hill *
-Hills::get(int n) {
+Hills::get(int n) const {
if (n < 0 || n >= num) {
return NULL;
} else {
diff --git a/src/ImageMetaData.H b/src/ImageMetaData.H
index 5796964..da63db6 100644
--- a/src/ImageMetaData.H
+++ b/src/ImageMetaData.H
@@ -15,6 +15,10 @@ class ImageMetaData {
double direction;
double nick;
double tilt;
+ double k0;
+ double k1;
+ double u0;
+ double v0;
double focal_length_35mm;
double scale;
int projection_type;
@@ -37,6 +41,7 @@ class ImageMetaData {
double get_tilt();
double get_focal_length_35mm();
int get_projection_type();
+ void get_distortion_params(double *_k0, double *_k1);
void set_longitude(double v);
void set_latitude(double v);
@@ -46,5 +51,6 @@ class ImageMetaData {
void set_tilt(double v);
void set_focal_length_35mm(double v);
void set_projection_type(int v);
+ void set_distortion_params(double _k0, double _k1);
};
#endif
diff --git a/src/ImageMetaData.cxx b/src/ImageMetaData.cxx
index 6519239..d59db9d 100644
--- a/src/ImageMetaData.cxx
+++ b/src/ImageMetaData.cxx
@@ -24,6 +24,10 @@ ImageMetaData::ImageMetaData() {
direction = 0.0;
nick = 0.0;
tilt = 0.0;
+ k0 = 0.0;
+ k1 = 0.0;
+ u0 = 0.0;
+ v0 = 0.0;
focal_length_35mm = 35.0;
scale = NAN;
projection_type = 0;
@@ -124,7 +128,7 @@ ImageMetaData::load_image_exif(char *name) {
#define GIPFEL_FORMAT_1 "gipfel: longitude %lf, latitude %lf, height %lf, direction %lf, nick %lf, tilt %lf, scale %lf, projection type %d"
-#define GIPFEL_FORMAT_2 "gipfel: longitude %lf, latitude %lf, height %lf, direction %lf, nick %lf, tilt %lf, focal_length_35mm %lf, projection type %d"
+#define GIPFEL_FORMAT_2 "gipfel: longitude %lf, latitude %lf, height %lf, direction %lf, nick %lf, tilt %lf, focal_length_35mm %lf, projection type %d, k0 %lf, k1 %lf"
int
ImageMetaData::load_image_jpgcom(char *name) {
@@ -133,9 +137,9 @@ ImageMetaData::load_image_jpgcom(char *name) {
pid_t pid;
int status;
char buf[1024];
- double lo, la, he, dir, ni, ti, fr;
+ double lo, la, he, dir, ni, ti, fr, _k0, _k1;
int pt;
- int ret = 1;
+ int n, ret = 1;
args[0] = "rdjpgcom";
args[1] = name;
@@ -145,8 +149,8 @@ ImageMetaData::load_image_jpgcom(char *name) {
if (p) {
while (fgets(buf, sizeof(buf), p) != NULL) {
- if (sscanf(buf, GIPFEL_FORMAT_2,
- &lo, &la, &he, &dir, &ni, &ti, &fr, &pt) >= 7) {
+ if ((n = sscanf(buf, GIPFEL_FORMAT_2,
+ &lo, &la, &he, &dir, &ni, &ti, &fr, &pt, &_k0, &_k1)) >= 8) {
longitude = lo;
latitude = la;
@@ -157,6 +161,11 @@ ImageMetaData::load_image_jpgcom(char *name) {
focal_length_35mm = fr;
projection_type = pt;
+ if (n >= 10) {
+ k0 = _k0;
+ k1 = _k1;
+ }
+
ret = 0;
break;
@@ -226,7 +235,8 @@ ImageMetaData::save_image_jpgcom(char *in_img, char *out_img) {
nick,
tilt,
focal_length_35mm,
- projection_type);
+ projection_type,
+ k0, k1, u0, v0);
// try to save gipfel data in JPEG comment section
args[0] = "wrjpgcom";
@@ -338,3 +348,15 @@ void
ImageMetaData::set_projection_type(int v) {
projection_type = v;
}
+
+void
+ImageMetaData::get_distortion_params(double *_k0, double *_k1) {
+ *_k0 = k0;
+ *_k1 = k1;
+}
+
+void
+ImageMetaData::set_distortion_params(double _k0, double _k1) {
+ k0 = _k0;
+ k1 = _k1;
+}
diff --git a/src/Makefile.am b/src/Makefile.am
index 8ba69c1..1cae9b7 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -7,6 +7,7 @@ gipfel_SOURCES = \
Panorama.cxx \
Projection.cxx \
ProjectionTangential.cxx \
+ ProjectionTangentialLSQ.cxx \
ProjectionSphaeric.cxx \
Hill.cxx \
Fl_Value_Dial.cxx \
@@ -24,6 +25,7 @@ noinst_HEADERS = \
Panorama.H \
Projection.H \
ProjectionTangential.H \
+ ProjectionTangentialLSQ.H \
ProjectionSphaeric.H \
Hill.H \
ViewParams.H \
diff --git a/src/Makefile.lsq_funcs b/src/Makefile.lsq_funcs
new file mode 100644
index 0000000..f117e1e
--- /dev/null
+++ b/src/Makefile.lsq_funcs
@@ -0,0 +1,4 @@
+all: lsq_funcs.c
+
+lsq_funcs.c: lsq_funcs.mac
+ maxima -b lsq_funcs.mac | grep "^static" > lsq_funcs.c
diff --git a/src/Panorama.H b/src/Panorama.H
index 81cdf8e..bfb39b8 100644
--- a/src/Panorama.H
+++ b/src/Panorama.H
@@ -23,7 +23,6 @@ class Panorama {
Hills *mountains;
Hills *close_mountains;
Hills *visible_mountains;
- ViewParams parms;
Projection::Projection *proj;
Projection::Projection_t projection_type;
@@ -62,6 +61,8 @@ class Panorama {
double pi_d, deg2rad;
public:
+ ViewParams parms;
+
Panorama();
~Panorama();
@@ -122,7 +123,7 @@ class Panorama {
double get_real_distance(Hill *m);
- int comp_params(Hill *m1, Hill *m2);
+ int comp_params(Hills *h);
int guess(Hills *p1, Hill *m1);
diff --git a/src/Panorama.cxx b/src/Panorama.cxx
index 4ecf2f4..503f373 100644
--- a/src/Panorama.cxx
+++ b/src/Panorama.cxx
@@ -11,6 +11,7 @@
#include "Panorama.H"
#include "ProjectionTangential.H"
+#include "ProjectionTangentialLSQ.H"
#include "ProjectionSphaeric.H"
#define EARTH_RADIUS 6371010.0
@@ -160,9 +161,10 @@ int
Panorama::guess(Hills *p, Hill *m1) {
Hill *p2, *m_tmp1, *m_tmp2;
Hill *m2;
+ Hills h;
double best = 100000000.0, v;
double a_center_best, a_nick_best, a_tilt_best, scale_best;
- int x1_sav, y1_sav;
+ double x1_sav, y1_sav;
int i, j;
if (m1 == NULL) {
@@ -192,7 +194,10 @@ Panorama::guess(Hills *p, Hill *m1) {
m2->x = p2->x;
m2->y = p2->y;
- comp_params(m1, m2);
+ h.clear();
+ h.add(m1);
+ h.add(m2);
+ comp_params(&h);
v = get_value(p);
@@ -219,10 +224,10 @@ Panorama::guess(Hills *p, Hill *m1) {
}
int
-Panorama::comp_params(Hill *m1, Hill *m2) {
+Panorama::comp_params(Hills *h) {
int ret;
- ret = proj->comp_params(m1, m2, &parms);
+ ret = proj->comp_params(h, &parms);
update_visible_mountains();
return ret;
}
@@ -285,7 +290,7 @@ Panorama::set_projection(Projection::Projection_t p) {
switch(projection_type) {
case Projection::TANGENTIAL:
- proj = new ProjectionTangential();
+ proj = new ProjectionTangentialLSQ();
view_angle = pi_d / 3.0;
break;
case Projection::SPHAERIC:
@@ -470,7 +475,6 @@ Panorama::update_visible_mountains() {
m = close_mountains->get(i);
m->a_view = m->alph - parms.a_center;
-
if (m->a_view > pi_d) {
m->a_view -= 2.0*pi_d;
} else if (m->a_view < -pi_d) {
@@ -494,11 +498,7 @@ Panorama::update_coordinates() {
for (int i=0; i<visible_mountains->get_num(); i++) {
m = visible_mountains->get(i);
- double tmp_x, tmp_y;
-
- proj->get_coordinates(m->a_view, m->a_nick, &parms, &tmp_x, &tmp_y);
- m->x = (int) rint(tmp_x);
- m->y = (int) rint(tmp_y);
+ proj->get_coordinates(m->alph, m->a_nick, &parms, &m->x, &m->y);
}
}
@@ -593,4 +593,3 @@ Panorama::get_coordinates(double a_view, double a_nick, double *x, double *y) {
return 0;
}
-
diff --git a/src/Projection.H b/src/Projection.H
index 3bc816e..2a953b2 100644
--- a/src/Projection.H
+++ b/src/Projection.H
@@ -25,6 +25,6 @@ class Projection {
virtual void get_coordinates(double a_view, double a_nick,
const ViewParams *parms, double *x, double *y);
- virtual int comp_params(const Hill *m1, const Hill *m2, ViewParams *parms);
+ virtual int comp_params(const Hills *h, ViewParams *parms);
};
#endif
diff --git a/src/Projection.cxx b/src/Projection.cxx
index 02674a7..bb577a0 100644
--- a/src/Projection.cxx
+++ b/src/Projection.cxx
@@ -20,7 +20,7 @@ Projection::get_coordinates(double a_view, double a_nick,
}
int
-Projection::comp_params(const Hill *m1, const Hill *m2, ViewParams *parms) {
+Projection::comp_params(const Hills *h, ViewParams *parms) {
fprintf(stderr, "Error: Projection::comp_params()\n");
return 1;
}
diff --git a/src/ProjectionSphaeric.H b/src/ProjectionSphaeric.H
index c6ce5f4..bb47fef 100644
--- a/src/ProjectionSphaeric.H
+++ b/src/ProjectionSphaeric.H
@@ -22,6 +22,6 @@ class ProjectionSphaeric : public Projection {
void get_coordinates(double a_view, double a_nick, const ViewParams *parms,
double *x, double *y);
- int comp_params(const Hill *m1, const Hill *m2, ViewParams *parms);
+ int comp_params(const Hills *h, ViewParams *parms);
};
#endif
diff --git a/src/ProjectionSphaeric.cxx b/src/ProjectionSphaeric.cxx
index 96af4e0..5227355 100644
--- a/src/ProjectionSphaeric.cxx
+++ b/src/ProjectionSphaeric.cxx
@@ -13,8 +13,8 @@
#define BEST_UNDEF 10000000.0
int
-ProjectionSphaeric::comp_params(const Hill *m1, const Hill *m2, ViewParams *parms) {
- const Hill *m_tmp;
+ProjectionSphaeric::comp_params(const Hills *h, ViewParams *parms) {
+ const Hill *m_tmp, *m1, *m2;
double tmp_x, tmp_y;
double val;
ViewParams best, tmp;
@@ -22,6 +22,14 @@ ProjectionSphaeric::comp_params(const Hill *m1, const Hill *m2, ViewParams *parm
double d_m1_2, d_m2_2, d_m1_m2_2;
int i, j;
+
+ if (h->get_num() != 2) {
+ return 1;
+ }
+
+ m1 = h->get(0);
+ m2 = h->get(1);
+
if (m1->x < m2->x) {
m_tmp = m1;
m1 = m2;
diff --git a/src/ProjectionTangential.H b/src/ProjectionTangential.H
index 13e354a..6734e88 100644
--- a/src/ProjectionTangential.H
+++ b/src/ProjectionTangential.H
@@ -24,6 +24,6 @@ class ProjectionTangential : public Projection {
void get_coordinates(double a_view, double a_nick, const ViewParams *parms,
double *x, double *y);
- int comp_params(const Hill *m1, const Hill *m2, ViewParams *parms);
+ int comp_params(const Hills *h, ViewParams *parms);
};
#endif
diff --git a/src/ProjectionTangential.cxx b/src/ProjectionTangential.cxx
index 8eed552..16488ad 100644
--- a/src/ProjectionTangential.cxx
+++ b/src/ProjectionTangential.cxx
@@ -29,10 +29,18 @@ comp_tilt(double tan_nick_view, double tan_dir_view, double n_scale,
double x, double y, double pi_d);
int
-ProjectionTangential::comp_params(const Hill *m1, const Hill *m2, ViewParams *parms) {
- const Hill *tmp;
+ProjectionTangential::comp_params(const Hills *h, ViewParams *parms) {
+ const Hill *tmp, *m1, *m2;
double a_center_tmp, scale_tmp, a_nick_tmp;
+ if (h->get_num() != 2) {
+ return 1;
+ }
+
+ m1 = h->get(0);
+ m2 = h->get(1);
+
+
scale_tmp = comp_scale(m1->alph, m2->alph, m1->x, m2->x);
if (isnan(scale_tmp) || scale_tmp < 100.0) {
// try again with mountains swapped
diff --git a/src/ProjectionTangentialLSQ.H b/src/ProjectionTangentialLSQ.H
new file mode 100644
index 0000000..dea10ea
--- /dev/null
+++ b/src/ProjectionTangentialLSQ.H
@@ -0,0 +1,31 @@
+//
+// Copyright 2006 Johannes Hofmann <Johannes.Hofmann@gmx.de>
+//
+// This software may be used and distributed according to the terms
+// of the GNU General Public License, incorporated herein by reference.
+
+#ifndef PROJECTIONTANGENTIALLSQ_H
+#define PROJECTIONTANGENTIALLSQ_H
+
+#include "Hill.H"
+#include "Projection.H"
+
+class ProjectionTangentialLSQ : public Projection {
+ private:
+ double comp_center_angle(double alph_a, double alph_b,
+ double d1, double d2);
+
+ double comp_scale(double alph_a, double alph_b, double d1, double d2);
+
+ double angle_dist(double a1, double a2);
+
+ int lsq(const Hills *m, ViewParams *parms);
+
+ public:
+ void get_coordinates(double a_view, double a_nick,
+ const ViewParams *parms, double *x, double *y);
+
+ int comp_params(const Hills *h, ViewParams *parms);
+
+};
+#endif
diff --git a/src/ProjectionTangentialLSQ.cxx b/src/ProjectionTangentialLSQ.cxx
new file mode 100644
index 0000000..9bbf81c
--- /dev/null
+++ b/src/ProjectionTangentialLSQ.cxx
@@ -0,0 +1,355 @@
+//
+// Copyright 2006 Johannes Hofmann <Johannes.Hofmann@gmx.de>
+//
+// This software may be used and distributed according to the terms
+// of the GNU General Public License, incorporated herein by reference.
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+
+#include <gsl/gsl_rng.h>
+#include <gsl/gsl_randist.h>
+#include <gsl/gsl_vector.h>
+#include <gsl/gsl_blas.h>
+#include <gsl/gsl_multifit_nlin.h>
+
+#include "ProjectionTangentialLSQ.H"
+
+
+static double sec(double a) {
+ return 1.0 / cos(a);
+}
+
+static double pi_d = asin(1.0) * 2.0, deg2rad = pi_d / 180.0;
+
+#include "lsq_funcs.c"
+
+static double
+comp_tilt(double tan_nick_view, double tan_dir_view, double n_scale,
+ double tan_nick_m, double tan_dir_m,
+ double x, double y, double pi_d);
+
+int
+ProjectionTangentialLSQ::comp_params(const Hills *h, ViewParams *parms) {
+ const Hill *tmp, *m1, *m2;
+ double a_center_tmp, scale_tmp, a_nick_tmp;
+
+ if (h->get_num() < 2) {
+ fprintf(stderr, "Please position at least 3 hills\n");
+ return 1;
+ } else if (h->get_num() > 3) {
+ fprintf(stderr, "Performing calibration\n");
+ parms->k0 = 0.0;
+ parms->k1 = 0.0;
+ }
+
+ m1 = h->get(0);
+ m2 = h->get(1);
+
+ scale_tmp = comp_scale(m1->alph, m2->alph, m1->x, m2->x);
+ if (isnan(scale_tmp) || scale_tmp < 100.0) {
+ // try again with mountains swapped
+ tmp = m1;
+ m1 = m2;
+ m2 = tmp;
+ scale_tmp = comp_scale(m1->alph, m2->alph, m1->x, m2->x);
+ }
+
+
+ if (isnan(scale_tmp) || scale_tmp < 100.0) {
+ return 1;
+ } else {
+
+ parms->a_center = m1->alph;
+ parms->scale = scale_tmp;
+ parms->a_nick = 0.0;
+ parms->a_tilt = 0.0;
+
+ lsq(h, parms);
+
+ return 0;
+ }
+}
+
+double
+ProjectionTangentialLSQ::angle_dist(double a1, double a2) {
+ double ret;
+
+ a1 = fmod(a1, 2.0 * pi_d);
+ if (a1 < 0.0) {
+ a1 = a1 + 2.0 * pi_d;
+ }
+ a2 = fmod(a2, 2.0 * pi_d);
+ if (a2 < 0.0) {
+ a2 = a2 + 2.0 * pi_d;
+ }
+
+ ret = fabs(a1 - a2);
+ if (ret > pi_d) {
+ ret = 2.0 * pi_d - ret;
+ }
+
+ return ret;
+}
+
+struct data {
+ const Hills *h;
+ const ViewParams *old_params;
+};
+
+#define CALL(A) A(c_view, c_nick, c_tilt, scale, k0, k1, m->alph, m->a_nick)
+
+static int
+lsq_f (const gsl_vector * x, void *data, gsl_vector * f) {
+ struct data *dat = (struct data *) data;
+ double c_view, c_nick, c_tilt, scale, k0, k1, u0, v0;
+
+ c_view = gsl_vector_get (x, 0);
+ c_nick = gsl_vector_get (x, 1);
+ c_tilt = gsl_vector_get (x, 2);
+ scale = gsl_vector_get (x, 3);
+ if (x->size >= 6) {
+ k0 = gsl_vector_get (x, 4);
+ k1 = gsl_vector_get (x, 5);
+ } else {
+ k0 = dat->old_params->k0;
+ k1 = dat->old_params->k1;
+ }
+
+ for (int i=0; i<dat->h->get_num(); i++) {
+ Hill *m = dat->h->get(i);
+
+ double mx = CALL(mac_x);
+ double my = CALL(mac_y);
+
+ gsl_vector_set (f, i*2, mx - m->x);
+ gsl_vector_set (f, i*2+1, my - m->y);
+ }
+
+ return GSL_SUCCESS;
+}
+
+
+static int
+lsq_df (const gsl_vector * x, void *data, gsl_matrix * J) {
+ struct data *dat = (struct data *) data;
+ double c_view, c_nick, c_tilt, scale, k0, k1, u0, v0;
+
+ c_view = gsl_vector_get (x, 0);
+ c_nick = gsl_vector_get (x, 1);
+ c_tilt = gsl_vector_get (x, 2);
+ scale = gsl_vector_get (x, 3);
+ if (x->size >= 6) {
+ k0 = gsl_vector_get (x, 4);
+ k1 = gsl_vector_get (x, 5);
+ } else {
+ k0 = dat->old_params->k0;
+ k1 = dat->old_params->k1;
+ }
+
+ for (int i=0; i<dat->h->get_num(); i++) {
+ Hill *m = dat->h->get(i);
+
+ gsl_matrix_set (J, 2*i, 0, CALL(mac_x_dc_view));
+ gsl_matrix_set (J, 2*i, 1, CALL(mac_x_dc_nick));
+ gsl_matrix_set (J, 2*i, 2, CALL(mac_x_dc_tilt));
+ gsl_matrix_set (J, 2*i, 3, CALL(mac_x_dscale));
+ if (x->size >= 6) {
+ gsl_matrix_set (J, 2*i, 4, CALL(mac_x_dk0));
+ gsl_matrix_set (J, 2*i, 5, CALL(mac_x_dk1));
+ }
+
+ gsl_matrix_set (J, 2*i+1, 0, CALL(mac_y_dc_view));
+ gsl_matrix_set (J, 2*i+1, 1, CALL(mac_y_dc_nick));
+ gsl_matrix_set (J, 2*i+1, 2, CALL(mac_y_dc_tilt));
+ gsl_matrix_set (J, 2*i+1, 3, CALL(mac_y_dscale));
+ if (x->size >= 6) {
+ gsl_matrix_set (J, 2*i+1, 4, CALL(mac_y_dk0));
+ gsl_matrix_set (J, 2*i+1, 5, CALL(mac_y_dk1));
+ }
+ }
+
+ return GSL_SUCCESS;
+}
+
+static int
+lsq_fdf (const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J) {
+ lsq_f (x, data, f);
+ lsq_df (x, data, J);
+
+ return GSL_SUCCESS;
+}
+
+int
+ProjectionTangentialLSQ::lsq(const Hills *h, ViewParams *parms) {
+ const gsl_multifit_fdfsolver_type *T;
+ gsl_multifit_fdfsolver *s;
+ gsl_multifit_function_fdf f;
+ struct data dat;
+ double x_init[8];
+ gsl_vector_view x;
+ int status;
+ int num_params = h->get_num()>3?6:4;
+
+ fprintf(stderr, "x %f, y %f\n",
+ h->get(0)->x,
+ h->get(0)->y);
+
+
+ dat.h = h;
+ dat.old_params = parms;
+
+ x_init[0] = parms->a_center;
+ x_init[1] = parms->a_nick;
+ x_init[2] = parms->a_tilt;
+ x_init[3] = parms->scale;
+ x_init[4] = parms->k0;
+ x_init[5] = parms->k1;
+
+ x = gsl_vector_view_array (x_init, num_params);
+
+ f.f = &lsq_f;
+ f.df = &lsq_df;
+ f.fdf = &lsq_fdf;
+ f.n = h->get_num() * 2;
+ f.p = num_params;
+ f.params = &dat;
+
+ T = gsl_multifit_fdfsolver_lmsder;
+ s = gsl_multifit_fdfsolver_alloc (T, h->get_num() * 2, num_params);
+ gsl_multifit_fdfsolver_set (s, &f, &x.vector);
+
+ for (int i=0; i<100; i++) {
+
+ status = gsl_multifit_fdfsolver_iterate (s);
+ if (status) {
+ fprintf(stderr, "gsl_multifit_fdfsolver_iterate: %d\n", status);
+ break;
+ }
+
+ fprintf(stderr, "%d, |f(x)| = %g\n", i, gsl_blas_dnrm2 (s->f));
+ }
+
+ parms->a_center = gsl_vector_get(s->x, 0);
+ parms->a_nick = gsl_vector_get(s->x, 1);
+ parms->a_tilt = gsl_vector_get(s->x, 2);
+ parms->scale = gsl_vector_get(s->x, 3);
+
+ if (num_params == 6) {
+ parms->k0 = gsl_vector_get(s->x, 4);
+ parms->k1 = gsl_vector_get(s->x, 5);
+ }
+
+ gsl_multifit_fdfsolver_free (s);
+
+ double t_x, t_y;
+ get_coordinates(h->get(0)->a_view, h->get(0)->a_nick, parms, &t_x, &t_y);
+ fprintf(stderr, "center %f, view %f, nick %f, x %f (%f), dx %f, y %f (%f), dy %f\n",
+ parms->a_center / deg2rad,
+ h->get(0)->a_view, h->get(0)->a_nick,
+ h->get(0)->x,
+ t_x,
+ h->get(0)->x - mac_x(parms->a_center,
+ parms->a_nick,
+ parms->a_tilt,
+ parms->scale,
+ parms->k0,
+ parms->k1,
+ h->get(0)->a_view, h->get(0)->a_nick),
+ h->get(0)->y,
+ t_y,
+ h->get(0)->y - mac_y(parms->a_center,
+ parms->a_nick,
+ parms->a_tilt,
+ parms->scale,
+ parms->k0,
+ parms->k1,
+ h->get(0)->a_view, h->get(0)->a_nick));
+
+
+
+ return 0;
+}
+
+void
+ProjectionTangentialLSQ::get_coordinates(double alph, double a_nick,
+ const ViewParams *parms, double *x, double *y) {
+
+ *x = mac_x(parms->a_center, parms->a_nick, parms->a_tilt, parms->scale,
+ parms->k0, parms->k1, alph, a_nick);
+ *y = mac_y(parms->a_center, parms->a_nick, parms->a_tilt, parms->scale,
+ parms->k0, parms->k1, alph, a_nick);
+}
+
+double
+ProjectionTangentialLSQ::comp_center_angle(double a1, double a2, double d1, double d2) {
+ double sign1 = 1.0;
+ double tan_acenter, tan_a1, tan_a2, a_center;
+
+ tan_a1 = tan(a1);
+ tan_a2 = tan(a2);
+
+ tan_acenter = (((pow(((pow((1.0 + (tan_a1 * tan_a2)), 2.0) * ((d1 * d1) + (d2 * d2))) + (2.0 * d1 * d2 * ((2.0 * ((tan_a2 * tan_a1) - (tan_a2 * tan_a2))) - ((tan_a1 * tan_a1) * (2.0 + (tan_a2 * tan_a2))) - 1.0))), (1.0 / 2.0)) * sign1) + ((1.0 - (tan_a1 * tan_a2)) * (d1 - d2))) / (2.0 * ((d2 * tan_a2) - (d1 * tan_a1))));
+
+ a_center = atan(tan_acenter);
+
+ if (a_center > 2.0 * pi_d) {
+ a_center = a_center - 2.0 * pi_d;
+ } else if (a_center < 0.0) {
+ a_center = a_center + 2.0 * pi_d;
+ }
+
+ // atan(tan_dir_view) is not the only possible solution.
+ // Choose the one which is close to m1->alph.
+ if (fabs(a_center - a1) > pi_d/2.0) {
+ a_center = a_center + pi_d;
+ }
+
+ return a_center;
+}
+
+double
+ProjectionTangentialLSQ::comp_scale(double a1, double a2, double d1, double d2) {
+ double sign1 = 1.0;
+ double sc, tan_a1, tan_a2;
+
+ tan_a1 = tan(a1);
+ tan_a2 = tan(a2);
+
+ sc = ((((1.0 + (tan_a1 * tan_a2)) * (d1 - d2)) - (sign1 * pow((((1.0 + pow((tan_a1 * tan_a2), 2.0)) * ((d1 * d1) + (d2 * d2))) + (2.0 * ((tan_a1 * tan_a2 * pow((d1 + d2), 2.0)) - (d1 * d2 * (((tan_a1 * tan_a1) * (2.0 + (tan_a2 * tan_a2))) + 1.0 + (2.0 * (tan_a2 * tan_a2))))))), (1.0 / 2.0)))) / (2.0 * (tan_a1 - tan_a2)));
+
+ return sc;
+}
+
+static double
+comp_tilt(double tan_nick_view, double tan_dir_view, double n_scale,
+ double tan_nick_m, double tan_dir_m,
+ double x, double y, double pi_d) {
+ double y_tmp, x_tmp, sin_a_tilt1, sin_a_tilt2, sin_a_tilt, res;
+
+ y_tmp = - (((tan_nick_view - tan_nick_m) * n_scale) /
+ (tan_nick_m * tan_nick_view + 1));
+ x_tmp = - (((tan_dir_view - tan_dir_m) * n_scale) /
+ (tan_dir_m * tan_dir_view + 1));
+
+
+ sin_a_tilt1 = - (y * - pow(x*x + y*y - y_tmp*y_tmp, 0.5) - x * y_tmp) /
+ (x*x + y*y);
+
+ sin_a_tilt2 = - (y * pow(x*x + y*y - y_tmp*y_tmp, 0.5) - x * y_tmp) /
+ (x*x + y*y);
+
+ sin_a_tilt = fabs(sin_a_tilt1) < fabs(sin_a_tilt2)?sin_a_tilt1:sin_a_tilt2;
+
+ res = asin(sin_a_tilt);
+
+ if (res > pi_d / 4.0) {
+ res = res - pi_d / 2.0;
+ } else if (res < -pi_d / 4.0) {
+ res = res + pi_d / 2.0;
+ }
+
+ return res;
+}
diff --git a/src/ViewParams.H b/src/ViewParams.H
index d321253..2048fe2 100644
--- a/src/ViewParams.H
+++ b/src/ViewParams.H
@@ -13,6 +13,8 @@ class ViewParams {
double scale;
double a_nick;
double a_tilt;
+ double k0;
+ double k1;
};
#endif
diff --git a/src/lsq_funcs.mac b/src/lsq_funcs.mac
new file mode 100644
index 0000000..de71198
--- /dev/null
+++ b/src/lsq_funcs.mac
@@ -0,0 +1,40 @@
+/*
+ * This is the basic pinhole projection model with distortion correction
+ */
+
+x : tan(m_view - c_view);
+y : tan(c_nick - m_nick);
+x_rot : y * sin(c_tilt) + x * cos(c_tilt);
+y_rot : y * cos(c_tilt) - x * sin(c_tilt);
+d : x_rot ^ 2 + y_rot ^ 2;
+dist_fact : d ^ 2 * k1 + d * k0;
+x_dist : x_rot * (1 + dist_fact) * scale;
+y_dist : y_rot * (1 + dist_fact) * scale;
+
+/*
+ * Some mangling for C code generation
+ */
+
+x_expand : trigexpand(x_dist);
+y_expand : trigexpand(y_dist);
+
+args: "(double c_view, double c_nick, double c_tilt, double scale, double k0, double k1, double m_view, double m_nick)";
+
+printfunc(name, expression) := sprint("static double", name, args, "{ return ", string(subst(pow, "^", expression)), ";}", "
+");
+
+printfunc("mac_x", x_expand);
+printfunc("mac_y", y_expand);
+
+printfunc("mac_x_dc_view", diff(x_expand, c_view));
+printfunc("mac_x_dc_nick", diff(x_expand, c_nick));
+printfunc("mac_x_dc_tilt", diff(x_expand, c_tilt));
+printfunc("mac_x_dscale", diff(x_expand, scale));
+printfunc("mac_x_dk0", diff(x_expand, k0));
+printfunc("mac_x_dk1", diff(x_expand, k1));
+printfunc("mac_y_dc_view", diff(y_expand, c_view));
+printfunc("mac_y_dc_nick", diff(y_expand, c_nick));
+printfunc("mac_y_dc_tilt", diff(y_expand, c_tilt));
+printfunc("mac_y_dscale", diff(y_expand, scale));
+printfunc("mac_y_dk0", diff(y_expand, k0));
+printfunc("mac_y_dk1", diff(y_expand, k1));
generated by cgit v1.2.3 (git 2.25.1) at 2024-08-21 22:51:21 +0000