Line data Source code
1 : /******************************************************************************
2 : * $Id$
3 : *
4 : * Project: libgeotiff
5 : * Purpose: Code to abstract translation between pixel/line and PCS
6 : * coordinates.
7 : * Author: Frank Warmerdam, warmerda@home.com
8 : *
9 : ******************************************************************************
10 : * Copyright (c) 1999, Frank Warmerdam
11 : *
12 : * Permission is hereby granted, free of charge, to any person obtaining a
13 : * copy of this software and associated documentation files (the "Software"),
14 : * to deal in the Software without restriction, including without limitation
15 : * the rights to use, copy, modify, merge, publish, distribute, sublicense,
16 : * and/or sell copies of the Software, and to permit persons to whom the
17 : * Software is furnished to do so, subject to the following conditions:
18 : *
19 : * The above copyright notice and this permission notice shall be included
20 : * in all copies or substantial portions of the Software.
21 : *
22 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
23 : * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
24 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
25 : * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
26 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
27 : * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
28 : * DEALINGS IN THE SOFTWARE.
29 : *****************************************************************************/
30 :
31 : #include <math.h>
32 : #include <stddef.h>
33 :
34 : #include "geotiff.h"
35 : #include "geo_tiffp.h" /* external TIFF interface */
36 : #include "geo_keyp.h" /* private interface */
37 : #include "geokeys.h"
38 :
39 : /************************************************************************/
40 : /* inv_geotransform() */
41 : /* */
42 : /* Invert a 6 term geotransform style matrix. */
43 : /************************************************************************/
44 :
45 0 : static int inv_geotransform( double *gt_in, double *gt_out )
46 :
47 : {
48 : /* we assume a 3rd row that is [0 0 1] */
49 :
50 : /* Compute determinate */
51 :
52 0 : const double det = gt_in[0] * gt_in[4] - gt_in[1] * gt_in[3];
53 :
54 0 : if( fabs(det) < 0.000000000000001 )
55 0 : return 0;
56 :
57 0 : const double inv_det = 1.0 / det;
58 :
59 : /* compute adjoint, and divide by determinate */
60 :
61 0 : gt_out[0] = gt_in[4] * inv_det;
62 0 : gt_out[3] = -gt_in[3] * inv_det;
63 :
64 0 : gt_out[1] = -gt_in[1] * inv_det;
65 0 : gt_out[4] = gt_in[0] * inv_det;
66 :
67 0 : gt_out[2] = ( gt_in[1] * gt_in[5] - gt_in[2] * gt_in[4]) * inv_det;
68 0 : gt_out[5] = (-gt_in[0] * gt_in[5] + gt_in[2] * gt_in[3]) * inv_det;
69 :
70 0 : return 1;
71 : }
72 :
73 : /************************************************************************/
74 : /* GTIFTiepointTranslate() */
75 : /************************************************************************/
76 :
77 : static
78 0 : int GTIFTiepointTranslate( int gcp_count, double * gcps_in, double * gcps_out,
79 : double x_in, double y_in,
80 : double *x_out, double *y_out )
81 :
82 : {
83 : (void) gcp_count;
84 : (void) gcps_in;
85 : (void) gcps_out;
86 : (void) x_in;
87 : (void) y_in;
88 : (void) x_out;
89 : (void) y_out;
90 :
91 : /* I would appreciate a _brief_ block of code for doing second order
92 : polynomial regression here! */
93 0 : return FALSE;
94 : }
95 :
96 :
97 : /************************************************************************/
98 : /* GTIFImageToPCS() */
99 : /************************************************************************/
100 :
101 : /**
102 : * Translate a pixel/line coordinate to projection coordinates.
103 : *
104 : * At this time this function does not support image to PCS translations for
105 : * tiepoints-only definitions, only pixelscale and transformation matrix
106 : * formulations.
107 : *
108 : * @param gtif The handle from GTIFNew() indicating the target file.
109 : * @param x A pointer to the double containing the pixel offset on input,
110 : * and into which the easting/longitude will be put on completion.
111 : * @param y A pointer to the double containing the line offset on input,
112 : * and into which the northing/latitude will be put on completion.
113 : *
114 : * @return TRUE if the transformation succeeds, or FALSE if it fails. It may
115 : * fail if the file doesn't have properly setup transformation information,
116 : * or it is in a form unsupported by this function.
117 : */
118 :
119 0 : int GTIFImageToPCS( GTIF *gtif, double *x, double *y )
120 :
121 : {
122 0 : tiff_t *tif=gtif->gt_tif;
123 :
124 : int tiepoint_count;
125 0 : double *tiepoints = 0;
126 0 : if (!(gtif->gt_methods.get)(tif, GTIFF_TIEPOINTS,
127 : &tiepoint_count, &tiepoints ))
128 0 : tiepoint_count = 0;
129 :
130 : int count;
131 0 : double *pixel_scale = 0;
132 0 : if (!(gtif->gt_methods.get)(tif, GTIFF_PIXELSCALE, &count, &pixel_scale ))
133 0 : count = 0;
134 :
135 : int transform_count;
136 0 : double *transform = 0;
137 0 : if (!(gtif->gt_methods.get)(tif, GTIFF_TRANSMATRIX,
138 : &transform_count, &transform ))
139 0 : transform_count = 0;
140 :
141 : /* -------------------------------------------------------------------- */
142 : /* If the pixelscale count is zero, but we have tiepoints use */
143 : /* the tiepoint based approach. */
144 : /* -------------------------------------------------------------------- */
145 0 : int res = FALSE;
146 0 : if( tiepoint_count > 6 && count == 0 )
147 : {
148 0 : res = GTIFTiepointTranslate( tiepoint_count / 6,
149 : tiepoints, tiepoints + 3,
150 : *x, *y, x, y );
151 : }
152 :
153 : /* -------------------------------------------------------------------- */
154 : /* If we have a transformation matrix, use it. */
155 : /* -------------------------------------------------------------------- */
156 0 : else if( transform_count == 16 )
157 : {
158 0 : const double x_in = *x;
159 0 : const double y_in = *y;
160 :
161 0 : *x = x_in * transform[0] + y_in * transform[1] + transform[3];
162 0 : *y = x_in * transform[4] + y_in * transform[5] + transform[7];
163 :
164 0 : res = TRUE;
165 : }
166 :
167 : /* -------------------------------------------------------------------- */
168 : /* For now we require one tie point, and a valid pixel scale. */
169 : /* -------------------------------------------------------------------- */
170 0 : else if( count < 3 || tiepoint_count < 6 )
171 : {
172 0 : res = FALSE;
173 : }
174 :
175 : else
176 : {
177 0 : *x = (*x - tiepoints[0]) * pixel_scale[0] + tiepoints[3];
178 0 : *y = (*y - tiepoints[1]) * (-1 * pixel_scale[1]) + tiepoints[4];
179 :
180 0 : res = TRUE;
181 : }
182 :
183 : /* -------------------------------------------------------------------- */
184 : /* Cleanup */
185 : /* -------------------------------------------------------------------- */
186 0 : if(tiepoints)
187 0 : _GTIFFree(tiepoints);
188 0 : if(pixel_scale)
189 0 : _GTIFFree(pixel_scale);
190 0 : if(transform)
191 0 : _GTIFFree(transform);
192 :
193 0 : return res;
194 : }
195 :
196 : /************************************************************************/
197 : /* GTIFPCSToImage() */
198 : /************************************************************************/
199 :
200 : /**
201 : * Translate a projection coordinate to pixel/line coordinates.
202 : *
203 : * At this time this function does not support PCS to image translations for
204 : * tiepoints-only based definitions, only matrix and pixelscale/tiepoints
205 : * formulations are supposed.
206 : *
207 : * @param gtif The handle from GTIFNew() indicating the target file.
208 : * @param x A pointer to the double containing the pixel offset on input,
209 : * and into which the easting/longitude will be put on completion.
210 : * @param y A pointer to the double containing the line offset on input,
211 : * and into which the northing/latitude will be put on completion.
212 : *
213 : * @return TRUE if the transformation succeeds, or FALSE if it fails. It may
214 : * fail if the file doesn't have properly setup transformation information,
215 : * or it is in a form unsupported by this function.
216 : */
217 :
218 0 : int GTIFPCSToImage( GTIF *gtif, double *x, double *y )
219 :
220 : {
221 0 : tiff_t *tif=gtif->gt_tif;
222 0 : int result = FALSE;
223 :
224 : /* -------------------------------------------------------------------- */
225 : /* Fetch tiepoints and pixel scale. */
226 : /* -------------------------------------------------------------------- */
227 0 : double *tiepoints = NULL;
228 : int tiepoint_count;
229 0 : if (!(gtif->gt_methods.get)(tif, GTIFF_TIEPOINTS,
230 : &tiepoint_count, &tiepoints ))
231 0 : tiepoint_count = 0;
232 :
233 : int count;
234 0 : double *pixel_scale = NULL;
235 0 : if (!(gtif->gt_methods.get)(tif, GTIFF_PIXELSCALE, &count, &pixel_scale ))
236 0 : count = 0;
237 :
238 0 : int transform_count = 0;
239 0 : double *transform = NULL;
240 0 : if (!(gtif->gt_methods.get)(tif, GTIFF_TRANSMATRIX,
241 : &transform_count, &transform ))
242 0 : transform_count = 0;
243 :
244 : /* -------------------------------------------------------------------- */
245 : /* If the pixelscale count is zero, but we have tiepoints use */
246 : /* the tiepoint based approach. */
247 : /* -------------------------------------------------------------------- */
248 0 : if( tiepoint_count > 6 && count == 0 )
249 : {
250 0 : result = GTIFTiepointTranslate( tiepoint_count / 6,
251 : tiepoints + 3, tiepoints,
252 : *x, *y, x, y );
253 : }
254 :
255 : /* -------------------------------------------------------------------- */
256 : /* Handle matrix - convert to "geotransform" format, invert and */
257 : /* apply. */
258 : /* -------------------------------------------------------------------- */
259 0 : else if( transform_count == 16 )
260 : {
261 0 : const double x_in = *x;
262 0 : const double y_in = *y;
263 : double gt_in[6];
264 0 : gt_in[0] = transform[0];
265 0 : gt_in[1] = transform[1];
266 0 : gt_in[2] = transform[3];
267 0 : gt_in[3] = transform[4];
268 0 : gt_in[4] = transform[5];
269 0 : gt_in[5] = transform[7];
270 :
271 : double gt_out[6];
272 0 : if( !inv_geotransform( gt_in, gt_out ) )
273 0 : result = FALSE;
274 : else
275 : {
276 0 : *x = x_in * gt_out[0] + y_in * gt_out[1] + gt_out[2];
277 0 : *y = x_in * gt_out[3] + y_in * gt_out[4] + gt_out[5];
278 :
279 0 : result = TRUE;
280 : }
281 : }
282 :
283 : /* -------------------------------------------------------------------- */
284 : /* For now we require one tie point, and a valid pixel scale. */
285 : /* -------------------------------------------------------------------- */
286 0 : else if( count >= 3 && tiepoint_count >= 6 )
287 : {
288 0 : *x = (*x - tiepoints[3]) / pixel_scale[0] + tiepoints[0];
289 0 : *y = (*y - tiepoints[4]) / (-1 * pixel_scale[1]) + tiepoints[1];
290 :
291 0 : result = TRUE;
292 : }
293 :
294 : /* -------------------------------------------------------------------- */
295 : /* Cleanup. */
296 : /* -------------------------------------------------------------------- */
297 0 : if(tiepoints)
298 0 : _GTIFFree(tiepoints);
299 0 : if(pixel_scale)
300 0 : _GTIFFree(pixel_scale);
301 0 : if(transform)
302 0 : _GTIFFree(transform);
303 :
304 0 : return result;
305 : }
|
