[GRASS-SVN] r44349 -
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify
svn_grass at osgeo.org
svn_grass at osgeo.org
Tue Nov 16 06:21:44 EST 2010
Author: mmetz
Date: 2010-11-16 03:21:44 -0800 (Tue, 16 Nov 2010)
New Revision: 44349
Added:
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/bilinear.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/bilinear_f.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cubic.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cubic_f.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/nearest.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/readcell.c
Removed:
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/matrix.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/perform.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ps_cp.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/rowcol.h
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/write.c
Modified:
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ask_files.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ask_wind.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/aver_z.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cp.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/defs.h
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/exec.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/get_wind.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/global.h
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/local_proto.h
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/main.c
grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/rectify.c
Log:
fix for mountainous areas, add resampling methods bilinear and bicubic
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ask_files.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ask_files.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ask_files.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -27,6 +27,15 @@
repeat = 0;
for (i = 0; i < NFILES; i++)
err[i] = "";
+ for (i = 0; i < group.group_ref.nfiles && i < NFILES; i++) {
+ if (strlen(group.group_ref.file[i].name) < 15)
+ sprintf(result[i], group.group_ref.file[i].name);
+ else
+ sprintf(result[i], "%14s", group.group_ref.file[i].name);
+ }
+ for (; i < NFILES; i++)
+ result[i][0] = 0;
+
f1 = 0;
for (f2 = f1; f1 < group.group_ref.nfiles; f1 = f2) {
any = 0;
@@ -34,9 +43,6 @@
V_clear();
V_line(0,
_("Please select the file(s) you wish to rectify by naming an output file"));
- if (!repeat)
- for (i = 0; i < NFILES; i++)
- result[i][0] = 0;
repeat = 0;
for (i = 0; f2 < group.group_ref.nfiles && i < NFILES; i++, f2++) {
name = group.group_ref.file[f2].name;
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ask_wind.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ask_wind.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ask_wind.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -3,7 +3,7 @@
#include <grass/gis.h>
#include <grass/vask.h>
-static int round(double *);
+static int rnd(double *);
static int visually_equal(double, double);
int ask_window(struct Cell_head *window)
@@ -13,14 +13,12 @@
short ok;
struct Cell_head minimal;
-
-
- round(&window->north);
- round(&window->south);
- round(&window->west);
- round(&window->east);
- round(&window->ew_res);
- round(&window->ns_res);
+ rnd(&window->north);
+ rnd(&window->south);
+ rnd(&window->west);
+ rnd(&window->east);
+ rnd(&window->ew_res);
+ rnd(&window->ns_res);
G_copy(&minimal, window, sizeof(minimal));
window->rows = 0;
@@ -88,12 +86,12 @@
if (!V_call())
exit(1);
- round(&window->north);
- round(&window->south);
- round(&window->east);
- round(&window->west);
- round(&window->ew_res);
- round(&window->ns_res);
+ rnd(&window->north);
+ rnd(&window->south);
+ rnd(&window->east);
+ rnd(&window->west);
+ rnd(&window->ew_res);
+ rnd(&window->ns_res);
if ((window->ns_res <= 0) || (window->ew_res <= 0)) {
fprintf(stderr, "Illegal resolution value(s)\n");
@@ -194,7 +192,7 @@
return strcmp(xs, ys) == 0;
}
-static int round(double *x)
+static int rnd(double *x)
{
char xs[40];
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/aver_z.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/aver_z.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/aver_z.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -3,6 +3,7 @@
int get_aver_elev(struct Ortho_Control_Points *cpz, double *aver_z)
{
double meanz;
+ double *cp = cpz->z2;
int i, n;
/* Need 1 control points */
@@ -18,7 +19,7 @@
continue;
n++;
- meanz += *((cpz->z2)++);
+ meanz += *(cp++);
G_debug(3, "In ortho meanz = %f", meanz);
}
@@ -26,12 +27,5 @@
G_debug(1, "In ortho aver_z = %f", *aver_z);
- /* reset pointers */
- for (i = 0; i < cpz->count; i++) {
- if (cpz->status[i] <= 0)
- continue;
- (cpz->z2)--;
- }
-
return 0;
}
Added: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/bilinear.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/bilinear.c (rev 0)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/bilinear.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -0,0 +1,59 @@
+/*
+ * Name
+ * bilinear.c -- use bilinear interpolation for given row, col
+ *
+ * Description
+ * bilinear interpolation for the given row, column indices.
+ * If the given row or column is outside the bounds of the input map,
+ * the point in the output map is set to NULL.
+ * If any of the 4 surrounding points to be used in the interpolation
+ * is NULL it is filled with is neighbor value
+ */
+
+#include <math.h>
+#include <grass/gis.h>
+#include "global.h"
+
+void p_bilinear(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index */
+ double *col_idx, /* column index */
+ struct Cell_head *cellhd /* information of output map */
+ )
+{
+ int row; /* row indices for interp */
+ int col; /* column indices for interp */
+ int i, j;
+ DCELL t, u; /* intermediate slope */
+ DCELL result; /* result of interpolation */
+ DCELL c[2][2];
+
+ /* cut indices to integer */
+ row = (int)floor(*row_idx - 0.5);
+ col = (int)floor(*col_idx - 0.5);
+
+ /* check for out of bounds - if out of bounds set NULL value and return */
+ if (row < 0 || row + 1 >= cellhd->rows || col < 0 || col + 1 >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ for (i = 0; i < 2; i++)
+ for (j = 0; j < 2; j++) {
+ const DCELL *cellp = CPTR(ibuffer, row + i, col + j);
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+ c[i][j] = *cellp;
+ }
+
+ /* do the interpolation */
+ t = *col_idx - 0.5 - col;
+ u = *row_idx - 0.5 - row;
+
+ result = G_interp_bilinear(t, u, c[0][0], c[0][1], c[1][0], c[1][1]);
+
+ G_set_raster_value_d(obufptr, result, cell_type);
+}
Added: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/bilinear_f.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/bilinear_f.c (rev 0)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/bilinear_f.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -0,0 +1,49 @@
+/*
+ * Name
+ * bilinear_f.c -- use bilinear interpolation with fallback for given row, col
+ *
+ * Description
+ * bilinear interpolation for the given row, column indices.
+ * If the interpolated value (but not the nearest) is null,
+ * it falls back to nearest neighbor.
+ */
+
+#include <grass/gis.h>
+#include <math.h>
+#include "global.h"
+
+void p_bilinear_f(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index */
+ double *col_idx, /* column index */
+ struct Cell_head *cellhd /* cell header of input layer */
+ )
+{
+ /* start nearest neighbor to do some basic tests */
+ int row, col; /* row/col of nearest neighbor */
+ DCELL *cellp, cell;
+
+ /* cut indices to integer */
+ row = (int)floor(*row_idx);
+ col = (int)floor(*col_idx);
+
+ /* check for out of bounds - if out of bounds set NULL value */
+ if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ cellp = CPTR(ibuffer, row, col);
+ /* if nearest is null, all the other interps will be null */
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+ cell = *cellp;
+
+ p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
+ /* fallback to nearest if bilinear is null */
+ if (G_is_d_null_value(obufptr))
+ G_set_raster_value_d(obufptr, cell, cell_type);
+}
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cp.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cp.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cp.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -12,24 +12,23 @@
/* compute photo coordinates of image control points */
/* I_convert_con_points (group.name, &group.control_points, &group.photo_points, group.E12, group.N12); */
-
- sprintf(msg, "Control Z Point file for group [%s] in [%s] \n \n",
+ sprintf(msg, _("Control Z Point file for group [%s] in [%s] \n \n"),
group.name, G_mapset());
- fprintf(stderr, "Computing equations...\n\n");
+ G_verbose_message(_("Computing equations..."));
Compute_ortho_equation();
switch (group.con_equation_stat) {
case -1:
- strcat(msg, "Poorly placed Control Points!\n");
- strcat(msg, "Can not generate the transformation equation.\n");
- strcat(msg, "Run OPTION 7 again!\n");
+ strcat(msg, _("Poorly placed Control Points!\n"));
+ strcat(msg, _("Can not generate the transformation equation.\n"));
+ strcat(msg, _("Run OPTION 7 again!\n"));
break;
case 0:
- strcat(msg, "No active Control Points!\n");
- strcat(msg, "Can not generate the transformation equation.\n");
- strcat(msg, "Run OPTION 7 !\n");
+ strcat(msg, _("No active Control Points!\n"));
+ strcat(msg, _("Can not generate the transformation equation.\n"));
+ strcat(msg, _("Run OPTION 7 !\n"));
break;
default:
return 1;
@@ -46,21 +45,21 @@
if (!I_get_ref_points(group.name, &group.photo_points))
exit(0);
- sprintf(msg, "Reference Point file for group [%s] in [%s] \n \n",
+ sprintf(msg, _("Reference Point file for group [%s] in [%s] \n \n"),
group.name, G_mapset());
Compute_ref_equation();
switch (group.ref_equation_stat) {
case -1:
- strcat(msg, "Poorly placed Reference Points!\n");
- strcat(msg, "Can not generate the transformation equation.\n");
- strcat(msg, "Run OPTION 5 again!\n");
+ strcat(msg, _("Poorly placed Reference Points!\n"));
+ strcat(msg, _("Can not generate the transformation equation.\n"));
+ strcat(msg, _("Run OPTION 5 again!\n"));
break;
case 0:
- strcat(msg, "No active Reference Points!\n");
- strcat(msg, "Can not generate the transformation equation.\n");
- strcat(msg, "Run OPTION 5!\n");
+ strcat(msg, _("No active Reference Points!\n"));
+ strcat(msg, _("Can not generate the transformation equation.\n"));
+ strcat(msg, _("Run OPTION 5!\n"));
break;
default:
E12a = E12[0];
Added: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cubic.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cubic.c (rev 0)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cubic.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -0,0 +1,67 @@
+/*
+ * Name
+ * cubic.c -- use cubic convolution interpolation for given row, col
+ *
+ * Description
+ * cubic returns the value in the buffer that is the result of cubic
+ * convolution interpolation for the given row, column indices.
+ * If the given row or column is outside the bounds of the input map,
+ * the corresponding point in the output map is set to NULL.
+ *
+ * If single surrounding points in the interpolation matrix are
+ * NULL they where filled with their neighbor
+ */
+
+#include <grass/gis.h>
+#include <math.h>
+#include "global.h"
+
+void p_cubic(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index (decimal) */
+ double *col_idx, /* column index (decimal) */
+ struct Cell_head *cellhd /* information of output map */
+ )
+{
+ int row; /* row indices for interp */
+ int col; /* column indices for interp */
+ int i, j;
+ DCELL t, u; /* intermediate slope */
+ DCELL result; /* result of interpolation */
+ DCELL val[4]; /* buffer for temporary values */
+ DCELL cell[4][4];
+
+ /* cut indices to integer */
+ row = (int)floor(*row_idx - 0.5);
+ col = (int)floor(*col_idx - 0.5);
+
+ /* check for out of bounds of map - if out of bounds set NULL value */
+ if (row - 1 < 0 || row + 2 >= cellhd->rows ||
+ col - 1 < 0 || col + 2 >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ for (i = 0; i < 4; i++)
+ for (j = 0; j < 4; j++) {
+ const DCELL *cellp = CPTR(ibuffer, row - 1 + i, col - 1 + j);
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+ cell[i][j] = *cellp;
+ }
+
+ /* do the interpolation */
+ t = *col_idx - 0.5 - col;
+ u = *row_idx - 0.5 - row;
+
+ for (i = 0; i < 4; i++) {
+ val[i] = G_interp_cubic(t, cell[i][0], cell[i][1], cell[i][2], cell[i][3]);
+ }
+
+ result = G_interp_cubic(u, val[0], val[1], val[2], val[3]);
+
+ G_set_raster_value_d(obufptr, result, cell_type);
+}
Added: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cubic_f.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cubic_f.c (rev 0)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/cubic_f.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -0,0 +1,54 @@
+/*
+ * Name
+ * cubic_f.c -- use cubic interpolation with fallback for given row, col
+ *
+ * Description
+ * cubic interpolation for the given row, column indices.
+ * If the interpolated value (but not the nearest) is null,
+ * it falls back to bilinear. If that interp value is null,
+ * it falls back to nearest neighbor.
+ */
+
+#include <grass/gis.h>
+#include <math.h>
+#include "global.h"
+
+void p_cubic_f(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index */
+ double *col_idx, /* column index */
+ struct Cell_head *cellhd /* cell header of input layer */
+ )
+{
+ /* start nearest neighbor to do some basic tests */
+ int row, col; /* row/col of nearest neighbor */
+ DCELL *cellp, cell;
+
+ /* cut indices to integer */
+ row = (int)floor(*row_idx);
+ col = (int)floor(*col_idx);
+
+ /* check for out of bounds - if out of bounds set NULL value */
+ if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ cellp = CPTR(ibuffer, row, col);
+ /* if nearest is null, all the other interps will be null */
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+ cell = *cellp;
+
+ p_cubic(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
+ /* fallback to bilinear if cubic is null */
+ if (G_is_d_null_value(obufptr)) {
+ p_bilinear(ibuffer, obufptr, cell_type, row_idx, col_idx, cellhd);
+ /* fallback to nearest if bilinear is null */
+ if (G_is_d_null_value(obufptr))
+ G_set_raster_value_d(obufptr, cell, cell_type);
+ }
+}
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/defs.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/defs.h 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/defs.h 2010-11-16 11:21:44 UTC (rev 44349)
@@ -1,7 +1,6 @@
#include <curses.h>
-#include "orthophoto.h"
#include <grass/rowio.h>
-#include "rowcol.h"
+#include "orthophoto.h"
/* this is a curses structure */
typedef struct
@@ -37,14 +36,6 @@
typedef struct
{
- double XT, YT, ZT; /* object space */
- int rowT, colT;
- double xt, yt; /* image space */
- int rowt, colt;
-} Tie_Point;
-
-typedef struct
-{
double E12[3], N12[3], E21[3], N21[3];
} Patch;
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/exec.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/exec.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/exec.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -23,16 +23,33 @@
struct History hist;
int colr_ok, hist_ok, cats_ok;
long start_time, rectify_time, compress_time;
+ double aver_z;
+ int elevfd;
+ struct cache *ebuffer;
+ G_debug(1, "Open elevation raster: ");
- /* allocate the output cell matrix */
- cell_buf = (CELL **) G_calloc(NROWS, sizeof(void *));
- n = NCOLS * G_raster_size(map_type);
- for (i = 0; i < NROWS; i++) {
- cell_buf[i] = (void *)G_malloc(n);
- G_set_null_value(cell_buf[i], NCOLS, map_type);
+ /* open elevation raster */
+ select_target_env();
+ G_set_window(&target_window);
+ G_debug(1, "target window: rs=%d cs=%d n=%f s=%f w=%f e=%f\n",
+ target_window.rows, target_window.cols, target_window.north,
+ target_window.south, target_window.west, target_window.east);
+
+ elevfd = G_open_cell_old(elev_layer, mapset_elev);
+ if (elevfd < 0) {
+ G_fatal_error(_("Could not open elevation raster"));
+ return 1;
}
+ G_debug(1, "elev layer = %s mapset elev = %s elevfd = %d",
+ elev_layer, mapset_elev, elevfd);
+ ebuffer = readcell(elevfd, seg_mb_elev, 1);
+ G_close_cell(elevfd);
+ /* get an average elevation of the control points */
+ /* this is used only if target cells have no elevation */
+ get_aver_elev(&group.control_points, &aver_z);
+
/* rectify each file */
for (n = 0; n < group.group_ref.nfiles; n++) {
G_debug(2, "I look for files to ortho rectify");
@@ -60,7 +77,7 @@
G_debug(2, "Starting the rectification...");
- if (rectify(name, mapset, result)) {
+ if (rectify(name, mapset, ebuffer, aver_z, result)) {
G_debug(2, "Done. Writing results...");
select_target_env();
if (cats_ok) {
@@ -75,10 +92,7 @@
G_write_history(result, &hist);
select_current_env();
time(&rectify_time);
- if (compress(result))
- time(&compress_time);
- else
- compress_time = rectify_time;
+ compress_time = rectify_time;
report(name, mapset, result, rectify_time - start_time,
compress_time - rectify_time, 1);
}
@@ -87,6 +101,8 @@
report(name, mapset, result, (long)0, (long)0, 0);
}
}
+ close(ebuffer->fd);
+ G_free(ebuffer);
G_done_msg(" ");
return 0;
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/get_wind.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/get_wind.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/get_wind.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -1,5 +1,6 @@
#include <stdlib.h>
#include <string.h>
+#include <math.h>
#include "global.h"
int get_target_window(void)
@@ -58,10 +59,12 @@
int georef_window(struct Cell_head *w1, struct Cell_head *w2)
{
- double n, e, z1;
+ double n, e, z1, ad;
double n0, e0;
double aver_z;
- double diffew, diffns;
+ struct _corner {
+ double n, e;
+ } nw, ne, se, sw;
/* get an average elevation from the active control points */
get_aver_elev(&group.control_points, &aver_z);
@@ -85,6 +88,8 @@
w2->north = w2->south = n;
w2->west = w2->east = e;
+ nw.n = n;
+ nw.e = e;
I_georef(w1->east, w1->north, &e0, &n0, group.E12, group.N12);
I_inverse_ortho_ref(e0, n0, aver_z, &e, &n, &z1, &group.camera_ref,
@@ -97,6 +102,8 @@
G_debug(1, "target x = %f y = %f", e, n);
+ ne.n = n;
+ ne.e = e;
if (n > w2->north)
w2->north = n;
if (n < w2->south)
@@ -116,6 +123,8 @@
w1->south, e0, n0);
G_debug(1, "target x = %f y = %f", e, n);
+ sw.n = n;
+ sw.e = e;
if (n > w2->north)
w2->north = n;
if (n < w2->south)
@@ -135,6 +144,8 @@
w1->south, e0, n0);
G_debug(1, "target x = %f y = %f", e, n);
+ se.n = n;
+ se.e = e;
if (n > w2->north)
w2->north = n;
if (n < w2->south)
@@ -144,20 +155,58 @@
if (e < w2->west)
w2->west = e;
+ /* this results in ugly res values, and ns_res != ew_res */
+ /* and is no good for rotation */
+ /*
w2->ns_res = (w2->north - w2->south) / w1->rows;
w2->ew_res = (w2->east - w2->west) / w1->cols;
+ */
- /* Miori Luca & Mauro Martinelli, ITC-irst 2003: extend region to
- * avoid cut-off of image edges in mountainous terrain:
- * extend target area by (empirically) 15%
- */
- diffew = (w2->east - w2->west);
- diffns = (w2->north - w2->south);
- w2->east = w2->east + 0.15 * diffew;
- w2->west = w2->west - 0.15 * diffew;
- w2->south = w2->south - 0.15 * diffns;
- w2->north = w2->north + 0.15 * diffns;
+ /* alternative: account for rotation and order > 1 */
+ /* N-S extends along western and eastern edge */
+ w2->ns_res = (sqrt((nw.n - sw.n) * (nw.n - sw.n) +
+ (nw.e - sw.e) * (nw.e - sw.e)) +
+ sqrt((ne.n - se.n) * (ne.n - se.n) +
+ (ne.e - se.e) * (ne.e - se.e))) / (2.0 * w1->rows);
+
+ /* E-W extends along northern and southern edge */
+ w2->ew_res = (sqrt((nw.n - ne.n) * (nw.n - ne.n) +
+ (nw.e - ne.e) * (nw.e - ne.e)) +
+ sqrt((sw.n - se.n) * (sw.n - se.n) +
+ (sw.e - se.e) * (sw.e - se.e))) / (2.0 * w1->cols);
+
+ /* make ew_res = ns_res */
+ w2->ns_res = (w2->ns_res + w2->ew_res) / 2.0;
+ w2->ew_res = w2->ns_res;
+
+ /* nice round values */
+ if (w2->ns_res > 1) {
+ if (w2->ns_res < 10) {
+ /* round to first decimal */
+ w2->ns_res = (int)(w2->ns_res * 10 + 0.5) / 10.0;
+ w2->ew_res = w2->ns_res;
+ }
+ else {
+ /* round to integer */
+ w2->ns_res = (int)(w2->ns_res + 0.5);
+ w2->ew_res = w2->ns_res;
+ }
+ }
+
+ /* adjust extends */
+ ad = w2->north > 0 ? 0.5 : -0.5;
+ w2->north = (int) (ceil(w2->north / w2->ns_res) + ad) * w2->ns_res;
+ ad = w2->south > 0 ? 0.5 : -0.5;
+ w2->south = (int) (floor(w2->south / w2->ns_res) + ad) * w2->ns_res;
+ ad = w2->east > 0 ? 0.5 : -0.5;
+ w2->east = (int) (ceil(w2->east / w2->ew_res) + ad) * w2->ew_res;
+ ad = w2->west > 0 ? 0.5 : -0.5;
+ w2->west = (int) (floor(w2->west / w2->ew_res) + ad) * w2->ew_res;
+
+ w2->rows = (w2->north - w2->south + w2->ns_res / 2.0) / w2->ns_res;
+ w2->cols = (w2->east - w2->west + w2->ew_res / 2.0) / w2->ew_res;
+
G_debug(1, "FINAL");
G_debug(1, "east = %f \n west = %f \n north = %f \n south = %f",
w2->east, w2->west, w2->north, w2->south);
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/global.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/global.h 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/global.h 2010-11-16 11:21:44 UTC (rev 44349)
@@ -7,43 +7,44 @@
#include <grass/imagery.h>
#include <grass/ortholib.h>
+#include <grass/glocale.h>
#include "defs.h"
-/* the large, the worse the results in mountaneous regions!! 128 is MAX!
- * but: the large, the slower - wants a dynamic implementation - TODO
- * possible solution: ratio local elevation range/camera height = 0.003
- */
-
-#define TIE_ROW_DIST 128
-#define TIE_COL_DIST 128
-
-#define NROWS 128
-#define NCOLS 128
-
-/* do not modify past this point */
-
#ifndef GLOBAL
#define GLOBAL extern
#endif
-#define IDX int
-
/* activate debug in Gmakefile */
#ifdef DEBUG3
GLOBAL FILE *Bugsr;
#endif
-GLOBAL ROWCOL row_map[NROWS][NCOLS];
-GLOBAL ROWCOL col_map[NROWS][NCOLS];
-GLOBAL ROWCOL row_min[NROWS];
-GLOBAL ROWCOL row_max[NROWS];
-GLOBAL ROWCOL row_left[NROWS];
-GLOBAL ROWCOL row_right[NROWS];
-GLOBAL IDX row_idx[NROWS];
-GLOBAL int matrix_rows, matrix_cols;
+#define L2BDIM 6
+#define BDIM (1<<(L2BDIM))
+#define L2BSIZE (2*(L2BDIM))
+#define BSIZE (1<<(L2BSIZE))
+#define HI(i) ((i)>>(L2BDIM))
+#define LO(i) ((i)&((BDIM)-1))
+typedef DCELL block[BDIM][BDIM];
+
+struct cache
+{
+ int fd;
+ int stride;
+ int nblocks;
+ block **grid;
+ block *blocks;
+ int *refs;
+};
+
+typedef void (*func) (struct cache *, void *, int, double *, double *, struct Cell_head *);
+
+GLOBAL func interpolate; /* interpolation routine */
+
+GLOBAL int seg_mb_img, seg_mb_elev;
+GLOBAL char *method;
GLOBAL int temp_fd;
-GLOBAL RASTER_MAP_TYPE map_type;
GLOBAL CELL **cell_buf;
GLOBAL char *temp_name;
@@ -58,8 +59,6 @@
GLOBAL struct Ortho_Camera_File_Ref cam_info;
GLOBAL struct Cell_head elevhd;
-GLOBAL DCELL *elevbuf;
-GLOBAL int elevfd;
GLOBAL char *elev_layer;
GLOBAL char *mapset_elev;
@@ -72,6 +71,4 @@
GLOBAL struct Cell_head target_window;
-GLOBAL Tie_Point **T_Point;
-
#include "local_proto.h"
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/local_proto.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/local_proto.h 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/local_proto.h 2010-11-16 11:21:44 UTC (rev 44349)
@@ -9,14 +9,14 @@
int ask_file_from_list(char *, char *);
/* ask_wind.c */
+int ask_method(void);
+
+/* ask_wind.c */
int ask_window(struct Cell_head *);
/* aver_z.c */
int get_aver_elev(struct Ortho_Control_Points *, double *);
-/* compress.c */
-int compress(char *);
-
/* conv.c */
int view_to_col(View *, int);
int view_to_row(View *, int);
@@ -50,24 +50,37 @@
int get_target_window(void);
int georef_window(struct Cell_head *, struct Cell_head *);
-/* matrix.c */
-int compute_georef_matrix(struct Cell_head *, struct Cell_head *);
+/* rectify.c */
+int rectify(char *, char *, struct cache *, double, char *);
-/* perform.c */
-int perform_georef(int, void *rast);
+/* readcell.c */
+struct cache *readcell(int, int, int);
+block *get_block(struct cache *, int);
-/* ps_cp.c */
-int get_psuedo_control_pt(int, int);
+#define BKIDX(c,y,x) ((y) * (c)->stride + (x))
+#define BKPTR(c,y,x) ((c)->grid[BKIDX((c),(y),(x))])
+#define BLOCK(c,y,x) (BKPTR((c),(y),(x)) ? BKPTR((c),(y),(x)) : get_block((c),BKIDX((c),(y),(x))))
+#define CPTR(c,y,x) (&(*BLOCK((c),HI((y)),HI((x))))[LO((y))][LO((x))])
-/* rectify.c */
-int rectify(char *, char *, char *);
-
/* report.c */
int report(char *, char *, char *, long, long, int);
/* target.c */
int get_target(char *);
-/* write.c */
-int write_map(char *);
-int write_matrix(int, int);
+/* declare resampling methods */
+/* bilinear.c */
+extern void p_bilinear(struct cache *, void *, int, double *, double *,
+ struct Cell_head *);
+/* cubic.c */
+extern void p_cubic(struct cache *, void *, int, double *, double *,
+ struct Cell_head *);
+/* nearest.c */
+extern void p_nearest(struct cache *, void *, int, double *, double *,
+ struct Cell_head *);
+/* bilinear_f.c */
+extern void p_bilinear_f(struct cache *, void *, int, double *, double *,
+ struct Cell_head *);
+/* cubic_f.c */
+extern void p_cubic_f(struct cache *, void *, int, double *, double *,
+ struct Cell_head *);
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/main.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/main.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/main.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -6,11 +6,12 @@
* Markus Neteler <neteler itc.it>,
* Bernhard Reiter <bernhard intevation.de>,
* Glynn Clements <glynn gclements.plus.com>,
- * Hamish Bowman <hamish_b yahoo.com>
+ * Hamish Bowman <hamish_b yahoo.com>,
+ * Markus Metz
*
* PURPOSE: Rectifies an image by using the image to photo coordinate
* transformation matrix
- * COPYRIGHT: (C) 1999-2008 by the GRASS Development Team
+ * COPYRIGHT: (C) 1999-2010 by the GRASS Development Team
*
* This program is free software under the GNU General Public
* License (>=v2). Read the file COPYING that comes with GRASS
@@ -66,7 +67,7 @@
elev_layer = (char *)G_malloc(GNAME_MAX * sizeof(char));
mapset_elev = (char *)G_malloc(GMAPSET_MAX * sizeof(char));
- /* get group ref */
+ /* find group */
strcpy(group.name, name);
if (!I_find_group(group.name))
G_fatal_error(_("Group [%s] not found"), group.name);
@@ -88,10 +89,9 @@
/* ask for files to be rectified */
ask_files(group.name);
-
G_debug(1, "Looking for elevation file in group: <%s>", group.name);
- /* get the block elevation layer raster map in target location */
+ /* get the block elevation layer raster map in target location */
if (!I_get_group_elev(group.name, elev_layer, mapset_elev, tl,
math_exp, units, nd))
G_fatal_error(_("No target elevation model selected for group <%s>"),
@@ -104,7 +104,7 @@
G_get_cellhd(elev_layer, mapset_elev, &elevhd);
select_current_env();
- /** look for camera info for this block **/
+ /** look for camera info for this block **/
if (!I_get_group_camera(group.name, camera))
G_fatal_error(_("No camera reference file selected for group <%s>"),
group.name);
@@ -116,7 +116,7 @@
/* get initial camera exposure station, if any */
if (I_find_initial(group.name)) {
if (!I_get_init_info(group.name, &group.camera_exp))
- G_warning(_("Bad format in initial exposusre station file for group <%s>"),
+ G_warning(_("Bad format in initial exposure station file for group <%s>"),
group.name);
}
@@ -130,14 +130,9 @@
select_current_env();
get_target_window();
- /* modify elevation values if needed */
+ /* ask for interpolation method and amount of memory to be used */
+ ask_method();
- /***
- select_target_env();
- ask_elev_data();
- select_current_env();
- ***/
-
/* go do it */
exec_rectify();
Deleted: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/matrix.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/matrix.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/matrix.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -1,97 +0,0 @@
-
-/**************************************************************
- * compute_georef_matrix (win1, win2)
- *
- */
-#include <stdlib.h>
-#include "global.h"
-
-static int cmp(const void *, const void *);
-
-int compute_georef_matrix(struct Cell_head *win1, struct Cell_head *win2)
-{
- ROWCOL *rmap, *cmap, rr, cc;
- int nrow1, nrow2;
- int ncol1, ncol2;
- double n1, w1, ns_res1, ew_res1;
- double n2, e2, ns_res2, ew_res2;
- double nx, ex;
- double NX, EX;
- int row, col;
- int min, max;
-
- ns_res1 = win1->ns_res;
- ew_res1 = win1->ew_res;
- nrow1 = win1->rows;
- ncol1 = win1->cols;
- n1 = win1->north;
- w1 = win1->west;
-
- ns_res2 = win2->ns_res;
- ew_res2 = win2->ew_res;
- nrow2 = win2->rows;
- ncol2 = win2->cols;
- n2 = win2->north;
- e2 = win2->west;
- matrix_rows = nrow2;
- matrix_cols = ncol2;
-
- /* georef equation is
- * ex = E21a + E21b * col + E21c * row
- * nx = N21a + N21b * col + N21c * row
- *
- * compute common code (for northing) outside east loop
- */
-
- for (n2 = win2->north, row = 0; row < nrow2; row++, n2 -= ns_res2) {
- rmap = row_map[row];
- cmap = col_map[row];
- min = max = -1;
-
-/* G_debug(3, "\n\t got row = \t%d", row); */
-
- /* compute common code */
- EX = E21a + E21c * row;
- NX = N21a + N21c * row;
-
- for (e2 = win2->west, col = 0; col < ncol2; col++, e2 += ew_res2) {
- /* georef e2,n2 */
- ex = EX + E21b * col;
- nx = NX + N21b * col;
-
- rr = (n1 - nx) / ns_res1;
- if (rr < 0 || rr >= nrow1)
- rr = -1;
- else if (min < 0)
- min = max = rr;
- else if (rr < min)
- min = rr;
- else if (rr > max)
- max = rr;
- *rmap++ = rr;
-
- cc = (ex - w1) / ew_res1;
- if (cc < 0 || cc >= ncol1)
- cc = -1;
- *cmap++ = cc;
-
-/* G_debug(4, "\n\tnx = \t%f \tex = \t%f \n\trr = \t%d \tcc = \t%d", nx, ex,rr, cc); */
- }
-
- row_min[row] = min;
- row_max[row] = max;
- row_left[row] = 0;
- row_right[row] = matrix_cols - 1;
- row_idx[row] = row;
- }
- qsort(row_idx, nrow2, sizeof(IDX), cmp);
-
- return 0;
-}
-
-static int cmp(const void *aa, const void *bb)
-{
- const IDX *a = aa, *b = bb;
-
- return (int)(row_min[*a] - row_min[*b]);
-}
Added: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/nearest.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/nearest.c (rev 0)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/nearest.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -0,0 +1,41 @@
+
+/*
+ * nearest.c - returns the nearest neighbor to a given
+ * x,y position
+ */
+
+#include <math.h>
+#include <grass/gis.h>
+#include "global.h"
+
+void p_nearest(struct cache *ibuffer, /* input buffer */
+ void *obufptr, /* ptr in output buffer */
+ int cell_type, /* raster map type of obufptr */
+ double *row_idx, /* row index in input matrix */
+ double *col_idx, /* column index in input matrix */
+ struct Cell_head *cellhd /* cell header of input layer */
+ )
+{
+ int row, col; /* row/col of nearest neighbor */
+ DCELL *cellp;
+
+ /* cut indices to integer and get nearest cell */
+ /* the row_idx, col_idx correction for bilinear/bicubic does not apply here */
+ row = (int)floor(*row_idx);
+ col = (int)floor(*col_idx);
+
+ /* check for out of bounds - if out of bounds set NULL value */
+ if (row < 0 || row >= cellhd->rows || col < 0 || col >= cellhd->cols) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ cellp = CPTR(ibuffer, row, col);
+
+ if (G_is_d_null_value(cellp)) {
+ G_set_null_value(obufptr, 1, cell_type);
+ return;
+ }
+
+ G_set_raster_value_d(obufptr, *cellp, cell_type);
+}
Deleted: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/perform.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/perform.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/perform.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -1,105 +0,0 @@
-
-/*======================================================================
- perform.c --
-
- The actual reading and writing of the source imagery
- into the target imagery.
-
-======================================================================*/
-
-#include "global.h"
-
-static ROWCOL *rmap, *cmap, left, right;
-static int do_cell(int, void *, void *);
-
-int rast_size;
-
-int perform_georef(int infd, void *rast)
-{
- int row;
- int curidx;
- int idx;
- int i;
-
- rast_size = G_raster_size(map_type);
- for (row = 0; row < matrix_rows; row++)
- G_set_null_value(cell_buf[row], matrix_cols, map_type);
-
- curidx = 0;
- while (1) {
- /* find first row */
- while (curidx < matrix_rows) {
- idx = row_idx[curidx];
- row = row_min[idx];
- if (row >= 0)
- break;
- curidx++;
- }
-
- /* G_debug(3, " curidx %d", curidx); */
-
- if (curidx >= matrix_rows)
- break;
-
- /* G_debug(3, "read row %d", row); */
-
- if (G_get_raster_row_nomask
- (infd, G_incr_void_ptr(rast, rast_size), row, map_type) < 0)
- return 0;
-
- for (i = curidx; i < matrix_rows; i++) {
- idx = row_idx[i];
- if (row != row_min[idx])
- break;
-
- /* G_debug(4, " process matrix row %d", idx); */
-
- rmap = row_map[idx];
- cmap = col_map[idx];
- left = row_left[idx];
- right = row_right[idx];
- do_cell(row, G_incr_void_ptr(rast, rast_size), cell_buf[idx]);
-
- row_min[idx]++;
- if (row_min[idx] > row_max[idx])
- row_min[idx] = -1;
- row_left[idx] = left;
- row_right[idx] = right;
- }
- }
-
- return 0;
-}
-
-static int do_cell(int row, void *in, void *out)
-{
- int col;
- void *inptr, *outptr;
-
- for (; left <= right; left++) {
- inptr = G_incr_void_ptr(in, cmap[left] * rast_size);
- outptr = G_incr_void_ptr(out, left * rast_size);
- if (rmap[left] < 0)
- continue;
- if (rmap[left] != row)
- break;
- G_raster_cpy(outptr, inptr, 1, map_type);
- }
- for (; left <= right; right--) {
- inptr = G_incr_void_ptr(in, cmap[right] * rast_size);
- outptr = G_incr_void_ptr(out, right * rast_size);
- if (rmap[right] < 0)
- continue;
- if (rmap[right] != row)
- break;
- G_raster_cpy(outptr, inptr, 1, map_type);
- }
- for (col = left; col <= right; col++) {
- inptr = G_incr_void_ptr(in, cmap[col] * rast_size);
- outptr = G_incr_void_ptr(out, col * rast_size);
- if (rmap[col] == row)
- G_raster_cpy(outptr, inptr, 1, map_type);
- }
-
- return 0;
-}
Deleted: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ps_cp.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ps_cp.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/ps_cp.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -1,82 +0,0 @@
-#include <stdlib.h>
-#include <string.h>
-#include <grass/glocale.h>
-#include "global.h"
-
-int get_psuedo_control_pt(int tie_row, int tie_col)
-{
- char msg[200];
- struct Ortho_Photo_Points ps_cp;
- int i, j, k;
-
- G_debug(1, "In ps_cp");
-
- /* allocate psuedo struct, max points are four */
- ps_cp.count = 4;
- ps_cp.e1 = (double *)G_malloc(4 * sizeof(double));
- ps_cp.n1 = (double *)G_malloc(4 * sizeof(double));
- ps_cp.e2 = (double *)G_malloc(4 * sizeof(double));
- ps_cp.n2 = (double *)G_malloc(4 * sizeof(double));
- ps_cp.status = (int *)G_malloc(4 * sizeof(int));
- G_debug(1, "ps_cp allocated");
-
- /* pseudo points are four corners taken from T_Points */
- k = 0;
- for (i = 0; i < 2; i++) {
- for (j = 0; j < 2; j++) {
- ps_cp.e1[k] = T_Point[tie_row + i][tie_col + j].xt;
- ps_cp.n1[k] = T_Point[tie_row + i][tie_col + j].yt;
- ps_cp.e2[k] = (j * ((T_Point[tie_row][tie_col + 1].XT -
- T_Point[tie_row][tie_col].XT) /
- target_window.ew_res));
- ps_cp.n2[k] =
- (i *
- ((T_Point[tie_row][tie_col].YT -
- T_Point[tie_row + 1][tie_col].YT) / target_window.ns_res));
- ps_cp.status[k] = 1;
-
- G_debug(2, "\t k = %d\t i = %d\t j = %d", k, i, j);
- G_debug(2, "\t\t e1[k] = %f", ps_cp.e1[k]);
- G_debug(2, "\t\t n1[k] = %f", ps_cp.n1[k]);
- G_debug(2, "\t\t e2[k] = %f", ps_cp.e2[k]);
- G_debug(2, "\t\t n2[k] = %f", ps_cp.n2[k]);
-
- k++;
- }
- }
-
- G_debug(1, "ps_cp initialized");
-
- switch (I_compute_ref_equations(&ps_cp, E12, N12, E21, N21)) {
- case -1:
- G_debug(1, "\tref_equ: case -1");
- strcat(msg, _("Poorly placed psuedo control points. "));
- strcat(msg, _("Cannot generate the transformation equation."));
- break;
- case 0:
- G_debug(1, "\tref_equ: case 0");
- strcat(msg, _("No active psuedo control points"));
- break;
- default:
- G_debug(1, "\tref equ: case good");
-
- E12a = E12[0];
- E12b = E12[1];
- E12c = E12[2];
- N12a = N12[0];
- N12b = N12[1];
- N12c = N12[2];
- E21a = E21[0];
- E21b = E21[1];
- E21c = E21[2];
- N21a = N21[0];
- N21b = N21[1];
- N21c = N21[2];
-
- G_debug(1, "\t\tE21 = %f\t %f\t %f", E21a, E21b, E21c);
- G_debug(1, "\t\tN21 = %f\t %f\t %f", N21a, N21b, N21c);
-
- return 1;
- }
- G_fatal_error(msg);
-}
Added: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/readcell.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/readcell.c (rev 0)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/readcell.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -0,0 +1,132 @@
+/*
+ * readcell.c - reads an entire cell layer into a buffer
+ *
+ */
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <sys/types.h>
+#include <sys/stat.h>
+#include <fcntl.h>
+#include <unistd.h>
+#include <grass/gis.h>
+#include <grass/glocale.h>
+#include "global.h"
+
+struct cache *readcell(int fdi, int size, int target_env)
+{
+ DCELL *tmpbuf;
+ struct cache *c;
+ int nrows;
+ int ncols;
+ int row;
+ char *filename;
+ int nx, ny;
+ int nblocks;
+ int i;
+
+ nrows = G_window_rows();
+ ncols = G_window_cols();
+
+ ny = (nrows + BDIM - 1) / BDIM;
+ nx = (ncols + BDIM - 1) / BDIM;
+
+ if (size > 0)
+ nblocks = size * ((1 << 20) / sizeof(block));
+ else
+ nblocks = (nx + ny) * 2; /* guess */
+
+ if (nblocks > nx * ny)
+ nblocks = nx * ny;
+
+ c = G_malloc(sizeof(struct cache));
+ c->stride = nx;
+ c->nblocks = nblocks;
+ c->grid = (block **) G_calloc(nx * ny, sizeof(block *));
+ c->blocks = (block *) G_malloc(nblocks * sizeof(block));
+ c->refs = (int *)G_malloc(nblocks * sizeof(int));
+
+ if (nblocks < nx * ny) {
+ /* Temporary file must be created in output location */
+ select_target_env();
+ filename = G_tempfile();
+ if (!target_env)
+ select_current_env();
+ c->fd = open(filename, O_RDWR | O_CREAT | O_EXCL, 0600);
+ if (c->fd < 0)
+ G_fatal_error(_("Unable to open temporary file"));
+ remove(filename);
+ }
+ else
+ c->fd = -1;
+
+ G_important_message(_("Allocating memory and reading input map..."));
+ G_percent(0, nrows, 5);
+
+ for (i = 0; i < c->nblocks; i++)
+ c->refs[i] = -1;
+
+ tmpbuf = (DCELL *) G_malloc(nx * sizeof(block));
+
+ for (row = 0; row < nrows; row += BDIM) {
+ int x, y;
+
+ for (y = 0; y < BDIM; y++) {
+ G_percent(row + y, nrows, 5);
+
+ if (row + y >= nrows)
+ break;
+
+ G_get_d_raster_row(fdi, &tmpbuf[y * nx * BDIM], row + y);
+ }
+
+ for (x = 0; x < nx; x++)
+ for (y = 0; y < BDIM; y++)
+ if (c->fd >= 0) {
+ if (write
+ (c->fd, &tmpbuf[(y * nx + x) * BDIM],
+ BDIM * sizeof(DCELL)) < 0)
+ G_fatal_error(_("Error writing segment file"));
+ }
+ else
+ memcpy(&c->blocks[BKIDX(c, HI(row), x)][LO(y)][0],
+ &tmpbuf[(y * nx + x) * BDIM],
+ BDIM * sizeof(DCELL));
+ }
+
+ G_free(tmpbuf);
+
+ if (c->fd < 0)
+ for (i = 0; i < c->nblocks; i++) {
+ c->grid[i] = &c->blocks[i];
+ c->refs[i] = i;
+ }
+
+ return c;
+}
+
+block *get_block(struct cache * c, int idx)
+{
+ int replace = rand() % c->nblocks;
+ block *p = &c->blocks[replace];
+ int ref = c->refs[replace];
+ off_t offset = (off_t) idx * sizeof(DCELL) << L2BSIZE;
+
+ if (c->fd < 0)
+ G_fatal_error(_("Internal error: cache miss on fully-cached map"));
+
+ if (ref >= 0)
+ c->grid[ref] = NULL;
+
+ c->grid[idx] = p;
+ c->refs[replace] = idx;
+
+ if (lseek(c->fd, offset, SEEK_SET) < 0)
+ G_fatal_error(_("Error seeking on segment file"));
+
+ if (read(c->fd, p, sizeof(block)) < 0)
+ G_fatal_error(_("Error writing segment file"));
+
+ return p;
+}
Modified: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/rectify.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/rectify.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/rectify.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -11,230 +11,116 @@
#include "global.h"
#include "local_proto.h"
-int rectify(char *name, char *mapset, char *result)
+int rectify(char *name, char *mapset, struct cache *ebuffer,
+ double aver_z, char *result)
{
- struct Cell_head cellhd, win;
+ struct Cell_head cellhd;
int ncols, nrows;
int row, col;
- int infd;
- void *rast;
- int x_ties, y_ties;
- int tie_row, tie_col;
- int i;
- double n2, e2, z2;
- double nx, ex, zx;
- int r2, c2;
- double row2, col2;
- double aver_z;
+ int infd, outfd;
+ RASTER_MAP_TYPE map_type;
+ int cell_size;
+ void *trast, *tptr;
+ double n1, e1, z1;
+ double nx, ex, nx1, ex1, zx1;
+ double row_idx, col_idx;
+ struct cache *ibuffer;
- G_debug(1, "Open temp elevation file: ");
+ select_current_env();
+ if (G_get_cellhd(name, mapset, &cellhd) < 0)
+ return 0;
- /* open temporary elevation cell layer */
- select_target_env();
+ /* open the file to be rectified
+ * set window to cellhd first to be able to read file exactly
+ */
+ G_set_window(&cellhd);
+ G_debug(2, "cellhd: rs=%d cs=%d n=%f s=%f w=%f e=%f\n",
+ cellhd.rows, cellhd.cols, cellhd.north,
+ cellhd.south, cellhd.west, cellhd.east);
- /**G_set_window (&elevhd);**/
- G_debug(1, "target window: rs=%d cs=%d n=%f s=%f w=%f e=%f\n",
- target_window.rows, target_window.cols, target_window.north,
- target_window.south, target_window.west, target_window.east);
-
- G_set_window(&target_window);
- elevfd = G_open_cell_old(elev_layer, mapset_elev);
-
- /**G_get_cellhd (elev_layer, mapset_elev, &elevhd);**/
- elevbuf = G_allocate_d_raster_buf(); /* enforce DCELL */
-
- /* get an average elevation of the control points */
- /* this is used only if TIE points are outside of the elev_layer boundary */
- get_aver_elev(&group.control_points, &aver_z);
-
-
- if (elevfd < 0) {
- G_debug(1, "CANT OPEN ELEV");
- G_debug(1, "elev layer = %s mapset elev = %s elevfd = %d",
- elev_layer, mapset_elev, elevfd);
+ infd = G_open_cell_old(name, mapset);
+ if (infd < 0) {
return 0;
}
+ map_type = G_get_raster_map_type(infd);
+ cell_size = G_raster_size(map_type);
+ if (strcmp(method, "nearest") != 0) {
+ map_type = DCELL_TYPE;
+ cell_size = G_raster_size(map_type);
+ }
- G_debug(1, "elev layer = %s mapset elev = %s elevfd = %d",
- elev_layer, mapset_elev, elevfd);
+ ibuffer = readcell(infd, seg_mb_img, 0);
- /* alloc Tie_Points */
- y_ties = (int)(target_window.rows / TIE_ROW_DIST) + 2;
- x_ties = (int)(target_window.cols / TIE_COL_DIST) + 2;
+ G_close_cell(infd); /* (pmx) 17 april 2000 */
- G_debug(1, "Number Tie_Points: y_ties %d \tx_ties %d", y_ties,
- x_ties);
+ G_message(_("Rectify <%s@%s> (location <%s>)"),
+ name, mapset, G_location());
+ select_target_env();
+ G_message(_("into <%s@%s> (location <%s>) ..."),
+ result, G_mapset(), G_location());
- T_Point = (Tie_Point **) G_malloc(y_ties * sizeof(Tie_Point *));
- for (i = 0; i < y_ties; i++)
- T_Point[i] = (Tie_Point *) G_malloc(x_ties * sizeof(Tie_Point));
+ G_set_window(&target_window);
+ nrows = target_window.rows;
+ ncols = target_window.cols;
- /* build Tie_Points */
- nrows = 0;
- for (tie_row = 0; tie_row < y_ties; tie_row++) {
- n2 = target_window.north -
- (tie_row * TIE_ROW_DIST * target_window.ns_res) - 1;
- if (n2 <= target_window.south)
- n2 = target_window.south + 1;
- row2 = northing_to_row(&target_window, n2);
- r2 = (int)row2;
+ outfd = G_open_raster_new(result, map_type);
+ trast = G_allocate_raster_buf(map_type);
- if ((G_get_d_raster_row(elevfd, elevbuf, r2)) < 0) {
- G_debug(1, "ERROR reading elevation layer %s fd = %d : row %d",
- elev_layer, elevfd, r2);
- exit(0);
- }
- ncols = 0;
- for (tie_col = 0; tie_col < x_ties; tie_col++) {
- e2 = target_window.west +
- (tie_col * TIE_COL_DIST * target_window.ew_res) + 1;
- if (e2 >= target_window.east)
- e2 = target_window.east - 1;
+ for (row = 0; row < nrows; row++) {
+ n1 = target_window.north - (row + 0.5) * target_window.ns_res;
- G_debug(5, "Tie_Point \t row %d \tcol %d", tie_row, tie_col);
- G_debug(5, "\t east %f\t north %f", e2, n2);
+ G_percent(row, nrows, 2);
- col2 = easting_to_col(&target_window, e2);
- c2 = (int)col2;
+ G_set_null_value(trast, ncols, map_type);
+ tptr = trast;
+ for (col = 0; col < ncols; col++) {
+ DCELL *zp = CPTR(ebuffer, row, col);
- G_debug(5, "\t\t row2 = %f \t col2 = %f", row2, col2);
- G_debug(5, "\t\t r2 = %d \t c2 = %d", r2, c2);
- G_debug(5, "\t\t elevbuf[c2] = %f", (DCELL) elevbuf[c2]);
-
- /* if target TIE point has no elevation, set to aver_z */
- if (G_is_d_null_value(&elevbuf[c2]))
- z2 = aver_z;
+ e1 = target_window.west + (col + 0.5) * target_window.ew_res;
+
+ /* if target cell has no elevation, set to aver_z */
+ if (G_is_d_null_value(zp)) {
+ G_warning(_("No elevation available at row = %d, col = %d"), row, col);
+ z1 = aver_z;
+ }
else
- z2 = (double)elevbuf[c2];
+ z1 = *zp;
- G_debug(5, "\t\t e2 = %f \t n2 = %f \t z2 = %f", e2, n2, z2);
- G_debug(5, "\t\t XC = %f \t YC = %f \t ZC = %f", group.XC,
- group.YC, group.ZC);
- G_debug(5, "\t\t omega = %f \t phi = %f \t kappa = %f",
- group.omega, group.phi, group.kappa);
-
- /* ex, nx: photo coordinates */
- I_ortho_ref(e2, n2, z2, &ex, &nx, &zx, &group.camera_ref,
+ /* target coordinates e1, n1 to photo coordinates ex1, nx1 */
+ I_ortho_ref(e1, n1, z1, &ex1, &nx1, &zx1, &group.camera_ref,
group.XC, group.YC, group.ZC, group.omega, group.phi,
group.kappa);
G_debug(5, "\t\tAfter ortho ref (photo cords): ex = %f \t nx = %f",
- ex, nx);
+ ex1, nx1);
- /* ex, nx: relative to (row,col) = 0 */
- I_georef(ex, nx, &ex, &nx, group.E21, group.N21);
+ /* photo coordinates ex1, nx1 to image coordinates ex, nx */
+ I_georef(ex1, nx1, &ex, &nx, group.E21, group.N21);
G_debug(5, "\t\tAfter geo ref: ex = %f \t nx = %f", ex, nx);
- T_Point[tie_row][tie_col].XT = e2;
- T_Point[tie_row][tie_col].YT = n2;
- T_Point[tie_row][tie_col].ZT = z2;
- T_Point[tie_row][tie_col].xt = ex;
- T_Point[tie_row][tie_col].yt = nx;
- }
+ /* convert to row/column indices of source raster */
+ row_idx = (cellhd.north - nx) / cellhd.ns_res;
+ col_idx = (ex - cellhd.west) / cellhd.ew_res;
- } /* end build */
+ /* resample data point */
+ interpolate(ibuffer, tptr, map_type, &row_idx, &col_idx, &cellhd);
- /* close elev layer so we can open the file to be rectified */
- select_target_env();
- if (!G_close_cell(elevfd)) {
- G_debug(1, "Can't close the elev file %s [%s in %s]",
- elev_layer, mapset_elev, G_location());
+ tptr = G_incr_void_ptr(tptr, cell_size);
+ }
+ G_put_raster_row(outfd, trast, map_type);
}
+ G_percent(1, 1, 1);
- /* open the result file into target window
- * this open must be first since we change the window later
- * raster maps open for writing are not affected by window changes
- * but those open for reading are
- *
- * also tell open that raster map will have the same format
- * (ie number of bytes per cell) as the file being rectified
- */
- select_current_env();
- if (G_get_cellhd(name, mapset, &cellhd) < 0)
- return 0;
+ G_free(trast);
- select_target_env();
- G_set_window(&target_window);
- G_set_cell_format(cellhd.format);
+ G_close_cell(outfd);
+ close(ibuffer->fd);
+ G_free(ibuffer);
+ if (G_get_cellhd(result, G_mapset(), &cellhd) < 0)
+ return 0;
- select_current_env();
-
- G_copy(&win, &target_window, sizeof(win));
-
- win.west += win.ew_res / 2;
- ncols = target_window.cols;
- col = 0;
-
- for (tie_col = 0; tie_col < (x_ties - 1); tie_col++) {
- G_debug(3, "Patching column %d:", ncols);
-
- if ((win.cols = ncols) > TIE_COL_DIST)
- win.cols = TIE_COL_DIST;
- win.north = target_window.north - win.ns_res / 2;
- nrows = target_window.rows;
- row = 0;
-
- for (tie_row = 0; tie_row < (y_ties - 1); tie_row++) {
- G_debug(5, "Patching %d row:", nrows);
-
- if ((win.rows = nrows) > TIE_ROW_DIST)
- win.rows = TIE_ROW_DIST;
-
- get_psuedo_control_pt(tie_row, tie_col);
-
- G_debug(5, "\t got psuedo pts: row %d \t col %d", tie_row, tie_col);
-
- compute_georef_matrix(&cellhd, &win);
-
- G_debug(5, "\t\tcompute geo matrix");
-
- /* open the source imagery file to be rectified */
- /* set window to cellhd first to be able to read file exactly */
- select_current_env();
- G_set_window(&cellhd);
- infd = G_open_cell_old(name, mapset);
- if (infd < 0) {
- close(infd);
- return 0;
- }
- map_type = G_get_raster_map_type(infd);
- rast =
- (void *)G_calloc(G_window_cols() + 1,
- G_raster_size(map_type));
- G_set_null_value(rast, G_window_cols() + 1, map_type);
-
- /* perform the actual data rectification */
- perform_georef(infd, rast);
-
- G_debug(5, "\t\tperform georef");
- G_debug(5, "\t\twrite matrix");
-
- /* close the source imagery file and free the buffer */
- select_current_env();
- G_close_cell(infd);
- G_free(rast);
- /* select_current_env(); */
- select_target_env();
-
-
- /* write of the data rectified into the result file */
- write_matrix(row, col);
-
- nrows -= win.rows;
- row += win.rows;
- win.north -= (win.ns_res * win.rows);
- }
-
- ncols -= win.cols;
- col += win.cols;
- win.west += (win.ew_res * win.cols);
- G_percent(col, col + ncols, 1);
- }
-
- select_target_env();
-
if (cellhd.proj == 0) { /* x,y imagery */
cellhd.proj = target_window.proj;
cellhd.zone = target_window.zone;
@@ -252,9 +138,6 @@
name, mapset);
}
- target_window.compressed = cellhd.compressed;
- G_close_cell(infd);
- write_map(result);
select_current_env();
return 1;
Deleted: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/rowcol.h
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/rowcol.h 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/rowcol.h 2010-11-16 11:21:44 UTC (rev 44349)
@@ -1 +0,0 @@
-#define ROWCOL short
Deleted: grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/write.c
===================================================================
--- grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/write.c 2010-11-16 11:16:22 UTC (rev 44348)
+++ grass/branches/develbranch_6/imagery/i.ortho.photo/photo.rectify/write.c 2010-11-16 11:21:44 UTC (rev 44349)
@@ -1,76 +0,0 @@
-#include <sys/types.h>
-#include <sys/stat.h>
-#include <fcntl.h>
-#include <unistd.h>
-#include <string.h>
-#include <errno.h>
-#include <grass/glocale.h>
-#include "global.h"
-
-
-int write_matrix(int row, int col)
-{
- int n;
-
- select_target_env();
-
- /* create temp file if it doesn't eexist */
- if (!temp_fd || (fcntl(temp_fd, F_GETFD) == -1)) {
- temp_name = G_tempfile();
- temp_fd = creat(temp_name, 0660);
- }
-
- for (n = 0; n < matrix_rows; n++) {
- off_t offset;
-
- offset =
- ((off_t) row++ * target_window.cols +
- col) * G_raster_size(map_type);
- lseek(temp_fd, offset, SEEK_SET);
-
- if (write(temp_fd, cell_buf[n], G_raster_size(map_type) * matrix_cols)
- != G_raster_size(map_type) * matrix_cols) {
- unlink(temp_name);
- G_fatal_error(_("Unable to write temp file: %s"),
- strerror(errno));
- }
- }
-
- select_current_env();
-
- return 0;
-}
-
-
-int write_map(char *name)
-{
- int fd, row;
- void *rast;
-
- G_set_window(&target_window);
-
- rast = G_allocate_raster_buf(map_type);
- close(temp_fd);
- temp_fd = open(temp_name, F_DUPFD);
- fd = G_open_raster_new(name, map_type);
-
- if (fd <= 0)
- G_fatal_error(_("Unable to open map %s"), name);
-
- for (row = 0; row < target_window.rows; row++) {
- if (read(temp_fd, rast, target_window.cols * G_raster_size(map_type))
- != target_window.cols * G_raster_size(map_type))
- G_fatal_error(_("Unable to write row %d"), row);
-
- if (G_put_raster_row(fd, rast, map_type) < 0) {
- G_fatal_error(_("Unable to write raster map. You might want to check available disk space and write permissions."));
- unlink(temp_name);
- }
- }
-
- close(temp_fd);
- unlink(temp_name);
- G_close_cell(fd);
-
- return 0;
-}
More information about the grass-commit
mailing list