1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: vtkGeoTransform.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15 /*-------------------------------------------------------------------------
16 Copyright 2008 Sandia Corporation.
17 Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
18 the U.S. Government retains certain rights in this software.
19 -------------------------------------------------------------------------*/
20
21 #include "vtkGeoTransform.h"
22
23 #include "vtkDoubleArray.h"
24 #include "vtkGeoProjection.h"
25 #include "vtkMath.h"
26 #include "vtkObjectFactory.h"
27 #include "vtkPoints.h"
28
29 #include "vtk_libproj.h"
30
31 vtkStandardNewMacro(vtkGeoTransform);
32 vtkCxxSetObjectMacro(vtkGeoTransform, SourceProjection, vtkGeoProjection);
33 vtkCxxSetObjectMacro(vtkGeoTransform, DestinationProjection, vtkGeoProjection);
34
vtkGeoTransform()35 vtkGeoTransform::vtkGeoTransform()
36 {
37 this->SourceProjection = nullptr;
38 this->DestinationProjection = nullptr;
39 }
40
~vtkGeoTransform()41 vtkGeoTransform::~vtkGeoTransform()
42 {
43 if ( this->SourceProjection )
44 {
45 this->SourceProjection->Delete();
46 }
47 if ( this->DestinationProjection )
48 {
49 this->DestinationProjection->Delete();
50 }
51 }
52
PrintSelf(ostream & os,vtkIndent indent)53 void vtkGeoTransform::PrintSelf( ostream& os, vtkIndent indent )
54 {
55 this->Superclass::PrintSelf( os, indent );
56 os << indent << "SourceProjection: " << this->SourceProjection << "\n";
57 os << indent << "DestinationProjection: " << this->DestinationProjection << "\n";
58 }
59
TransformPoints(vtkPoints * srcPts,vtkPoints * dstPts)60 void vtkGeoTransform::TransformPoints( vtkPoints* srcPts, vtkPoints* dstPts )
61 {
62 if ( ! srcPts || ! dstPts )
63 {
64 return;
65 }
66
67 vtkDoubleArray* srcCoords = vtkArrayDownCast<vtkDoubleArray>( srcPts->GetData() );
68 vtkDoubleArray* dstCoords = vtkArrayDownCast<vtkDoubleArray>( dstPts->GetData() );
69 if ( ! srcCoords || ! dstCoords )
70 { // data not in a form we can use directly anyway...
71 this->Superclass::TransformPoints( srcPts, dstPts );
72 return;
73 }
74 dstCoords->DeepCopy( srcCoords );
75
76 projPJ src = this->SourceProjection ? this->SourceProjection->GetProjection() : nullptr;
77 projPJ dst = this->DestinationProjection ? this->DestinationProjection->GetProjection() : nullptr;
78 if ( ! src && ! dst )
79 {
80 // we've already copied srcCoords to dstCoords and src=dst=0 implies no transform...
81 return;
82 }
83
84 if ( srcCoords->GetNumberOfComponents() < 2 )
85 {
86 vtkErrorMacro( << "Source coordinate array " << srcCoords << " only has " << srcCoords->GetNumberOfComponents()
87 << " components and at least 2 are required for geographic projections." );
88 return;
89 }
90
91 this->InternalTransformPoints( dstCoords->GetPointer( 0 ), dstCoords->GetNumberOfTuples(), dstCoords->GetNumberOfComponents() );
92 }
93
Inverse()94 void vtkGeoTransform::Inverse()
95 {
96 vtkGeoProjection* tmp = this->SourceProjection;
97 this->SourceProjection = this->DestinationProjection;
98 this->DestinationProjection = tmp;
99 this->Modified();
100 }
101
InternalTransformPoint(const float in[3],float out[3])102 void vtkGeoTransform::InternalTransformPoint( const float in[3], float out[3] )
103 {
104 double ind[3];
105 double oud[3];
106 int i;
107 for ( i = 0; i < 3; ++ i )
108 ind[i] = in[i];
109 this->InternalTransformPoint( ind, oud );
110 for ( i = 0; i < 3; ++ i )
111 out[i] = static_cast<float>(oud[i]);
112 }
113
InternalTransformPoint(const double in[3],double out[3])114 void vtkGeoTransform::InternalTransformPoint( const double in[3], double out[3] )
115 {
116 for ( int i = 0; i < 3; ++ i )
117 {
118 out[i] = in[i];
119 }
120 this->InternalTransformPoints( out, 1, 3 );
121 }
122
InternalTransformDerivative(const float in[3],float out[3],float derivative[3][3])123 void vtkGeoTransform::InternalTransformDerivative( const float in[3], float out[3], float derivative[3][3] )
124 {
125 double ind[3];
126 double oud[3];
127 double drd[3][3];
128 int i;
129 for ( i = 0; i < 3; ++ i )
130 ind[i] = in[i];
131 this->InternalTransformDerivative( ind, oud, drd );
132 for ( i = 0; i < 3; ++ i )
133 {
134 out[i] = static_cast<float>(oud[i]);
135 for ( int j = 0; j < 3; ++ j )
136 {
137 derivative[i][j] = drd[i][j];
138 }
139 }
140 }
141
InternalTransformDerivative(const double in[3],double out[3],double derivative[3][3])142 void vtkGeoTransform::InternalTransformDerivative( const double in[3], double out[3], double derivative[3][3] )
143 {
144 // FIXME: Need to use pj_factors for both source and inverted dest projection
145 (void) in;
146 (void) out;
147 (void) derivative;
148 }
149
150
MakeTransform()151 vtkAbstractTransform* vtkGeoTransform::MakeTransform()
152 {
153 vtkGeoTransform* geoTrans = vtkGeoTransform::New();
154 return geoTrans;
155 }
156
InternalTransformPoints(double * x,vtkIdType numPts,int stride)157 void vtkGeoTransform::InternalTransformPoints( double* x, vtkIdType numPts, int stride )
158 {
159 projPJ src = this->SourceProjection ? this->SourceProjection->GetProjection() : nullptr;
160 projPJ dst = this->DestinationProjection ? this->DestinationProjection->GetProjection() : nullptr;
161 int delta = stride - 2;
162 projLP lp;
163 projXY xy;
164 if ( src )
165 {
166 // Convert from src system to lat/long using inverse of src transform
167 double* coord = x;
168 for ( vtkIdType i = 0; i < numPts; ++ i )
169 {
170 xy.u = coord[0]; xy.v = coord[1];
171 lp = pj_inv( xy, src );
172 coord[0] = lp.u; coord[1] = lp.v;
173 coord += stride;
174 }
175 }
176 else // ! src
177 {
178 // src coords are in degrees, convert to radians
179 double* coord = x;
180 for ( vtkIdType i = 0; i < numPts; ++ i )
181 {
182 for ( int j = 0; j < 2; ++ j, ++ coord )
183 {
184 *coord = vtkMath::RadiansFromDegrees( *coord );
185 }
186 coord += delta;
187 }
188 }
189 if ( dst )
190 {
191 double* coord = x;
192 for ( vtkIdType i = 0; i < numPts; ++ i )
193 {
194 lp.u = coord[0]; lp.v = coord[1];
195 xy = pj_fwd( lp, dst );
196 coord[0] = xy.u; coord[1] = xy.v;
197 coord += stride;
198 }
199 }
200 else // ! dst
201 {
202 // dst coords are in radians, convert to degrees
203 double* coord = x;
204 for ( vtkIdType i = 0; i < numPts; ++ i )
205 {
206 for ( int j = 0; j < 2; ++ j, ++ coord )
207 {
208 *coord = vtkMath::DegreesFromRadians( *coord );
209 }
210 coord += delta;
211 }
212 }
213 }
214