Merge pull request #1322 from rouault/ossfuzz_fixes

Various ossfuzz fixes

11 files changed, 118 insertions, 23 deletions

diff --git a/src/ell_set.cpp b/src/ell_set.cpp
index 71746895..386b4f46 100644
--- a/src/ell_set.cpp
+++ b/src/ell_set.cpp
@@ -399,6 +399,10 @@ static int ellps_spherification (PJ *P) {
         break;
     }
 
+    if (P->a <= 0.) {
+        return proj_errno_set(P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
+    }
+
     /* Clean up the ellipsoidal parameters to reflect the sphere */
     P->es = P->e = P->f = 0;
     P->rf = HUGE_VAL;
diff --git a/src/gauss.cpp b/src/gauss.cpp
index b7a27191..49ccfa1c 100644
--- a/src/gauss.cpp
+++ b/src/gauss.cpp
@@ -65,9 +65,14 @@ void *pj_gauss_ini(double e, double phi0, double *chi, double *rc) {
     }
     *chi = asin(sphi / en->C);
     en->ratexp = 0.5 * en->C * e;
-    en->K = tan(.5 * *chi + M_FORTPI) / (
-        pow(tan(.5 * phi0 + M_FORTPI), en->C) *
-        srat(en->e * sphi, en->ratexp)  );
+    if( .5 * phi0 + M_FORTPI < 1e-10 ) {
+        en->K = 1.0 / srat(en->e * sphi, en->ratexp);
+    }
+    else {
+        en->K = tan(.5 * *chi + M_FORTPI) / (
+            pow(tan(.5 * phi0 + M_FORTPI), en->C) *
+            srat(en->e * sphi, en->ratexp)  );
+    }
     return ((void *)en);
 }
 
diff --git a/src/projections/aea.cpp b/src/projections/aea.cpp
index 9a0c4656..f457e836 100644
--- a/src/projections/aea.cpp
+++ b/src/projections/aea.cpp
@@ -155,6 +155,8 @@ static PJ *setup(PJ *P) {
     P->inv = e_inverse;
     P->fwd = e_forward;
 
+    if (fabs(Q->phi1) > M_HALFPI || fabs(Q->phi2) > M_HALFPI)
+        return destructor(P, PJD_ERR_LAT_LARGER_THAN_90);
     if (fabs(Q->phi1 + Q->phi2) < EPS10)
         return destructor(P, PJD_ERR_CONIC_LAT_EQUAL);
     Q->n = sinphi = sin(Q->phi1);
diff --git a/src/projections/hammer.cpp b/src/projections/hammer.cpp
index aa7d1ba9..b2e56a2c 100644
--- a/src/projections/hammer.cpp
+++ b/src/projections/hammer.cpp
@@ -24,7 +24,14 @@ static PJ_XY s_forward (PJ_LP lp, PJ *P) {           /* Spheroidal, forward */
     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
     double cosphi, d;
 
-    d = sqrt(2./(1. + (cosphi = cos(lp.phi)) * cos(lp.lam *= Q->w)));
+    cosphi = cos(lp.phi);
+    lp.lam *= Q->w;
+    double denom = 1. + cosphi * cos(lp.lam);
+    if( denom == 0.0 ) {
+        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
+        return proj_coord_error().xy;
+    }
+    d = sqrt(2./denom);
     xy.x = Q->m * d * cosphi * sin(lp.lam);
     xy.y = Q->rm * d * sin(lp.phi);
     return xy;
diff --git a/src/projections/ocea.cpp b/src/projections/ocea.cpp
index 75aa6666..4e28f727 100644
--- a/src/projections/ocea.cpp
+++ b/src/projections/ocea.cpp
@@ -15,8 +15,6 @@ struct pj_opaque {
     double  rtk;
     double  sinphi;
     double  cosphi;
-    double  singam;
-    double  cosgam;
 };
 } // anonymous namespace
 
@@ -61,15 +59,16 @@ PJ *PROJECTION(ocea) {
 
     Q->rok = 1. / P->k0;
     Q->rtk = P->k0;
+    double lam_p, phi_p;
     /*If the keyword "alpha" is found in the sentence then use 1point+1azimuth*/
     if ( pj_param(P->ctx, P->params, "talpha").i) {
         /*Define Pole of oblique transformation from 1 point & 1 azimuth*/
         alpha   = pj_param(P->ctx, P->params, "ralpha").f;
         lonz = pj_param(P->ctx, P->params, "rlonc").f;
         /*Equation 9-8 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/
-        Q->singam = atan(-cos(alpha)/(-sin(phi_0) * sin(alpha))) + lonz;
+        lam_p = atan(-cos(alpha)/(-sin(phi_0) * sin(alpha))) + lonz;
         /*Equation 9-7 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/
-        Q->sinphi = asin(cos(phi_0) * sin(alpha));
+        phi_p = asin(cos(phi_0) * sin(alpha));
     /*If the keyword "alpha" is NOT found in the sentence then use 2points*/
     } else {
         /*Define Pole of oblique transformation from 2 points*/
@@ -78,23 +77,30 @@ PJ *PROJECTION(ocea) {
         lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
         lam_2 = pj_param(P->ctx, P->params, "rlon_2").f;
         /*Equation 9-1 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/
-        Q->singam = atan2(cos(phi_1) * sin(phi_2) * cos(lam_1) -
+        lam_p = atan2(cos(phi_1) * sin(phi_2) * cos(lam_1) -
             sin(phi_1) * cos(phi_2) * cos(lam_2),
             sin(phi_1) * cos(phi_2) * sin(lam_2) -
             cos(phi_1) * sin(phi_2) * sin(lam_1) );
 
         /* take care of P->lam0 wrap-around when +lam_1=-90*/
         if (lam_1 == -M_HALFPI)
-            Q->singam = -Q->singam;
+            lam_p = -lam_p;
 
         /*Equation 9-2 page 80 (http://pubs.usgs.gov/pp/1395/report.pdf)*/
-        Q->sinphi = atan(-cos(Q->singam - lam_1) / tan(phi_1));
+        double cos_lamp_m_minus_lam_1 = cos(lam_p - lam_1);
+        double tan_phi_1 = tan(phi_1);
+        if( tan_phi_1 == 0.0 ) {
+            // Not sure if we want to support this case, but at least this avoids
+            // a division by zero, and gives the same result as the below atan()
+            phi_p = (cos_lamp_m_minus_lam_1 >= 0.0 ) ? -M_HALFPI : M_HALFPI;
+        }
+        else {
+            phi_p = atan(- cos_lamp_m_minus_lam_1 / tan_phi_1);
+        }
     }
-    P->lam0 = Q->singam + M_HALFPI;
-    Q->cosphi = cos(Q->sinphi);
-    Q->sinphi = sin(Q->sinphi);
-    Q->cosgam = cos(Q->singam);
-    Q->singam = sin(Q->singam);
+    P->lam0 = lam_p + M_HALFPI;
+    Q->cosphi = cos(phi_p);
+    Q->sinphi = sin(phi_p);
     P->inv = s_inverse;
     P->fwd = s_forward;
     P->es = 0.;
diff --git a/src/strerrno.cpp b/src/strerrno.cpp
index 5dedef98..a9310a55 100644
--- a/src/strerrno.cpp
+++ b/src/strerrno.cpp
@@ -30,7 +30,7 @@ pj_err_list[] = {
     "acos/asin: |arg| >1.+1e-14",                                      /* -19 */
     "tolerance condition error",                                       /* -20 */
     "conic lat_1 = -lat_2",                                            /* -21 */
-    "lat_1 >= 90",                                                     /* -22 */
+    "lat_0, lat_1 or lat_2 >= 90",                                     /* -22 */
     "lat_1 = 0",                                                       /* -23 */
     "lat_ts >= 90",                                                    /* -24 */
     "no distance between control points",                              /* -25 */
diff --git a/src/transformations/helmert.cpp b/src/transformations/helmert.cpp
index 034f76f4..c00fad2f 100644
--- a/src/transformations/helmert.cpp
+++ b/src/transformations/helmert.cpp
@@ -567,14 +567,14 @@ PJ *TRANSFORMATION(helmert, 0) {
     if (pj_param_exists (P->params, "theta")) {
         P->left  = PJ_IO_UNITS_PROJECTED;
         P->right = PJ_IO_UNITS_PROJECTED;
+        P->fwd    = helmert_forward;
+        P->inv    = helmert_reverse;
     }
 
     P->fwd4d  = helmert_forward_4d;
     P->inv4d  = helmert_reverse_4d;
     P->fwd3d  = helmert_forward_3d;
     P->inv3d  = helmert_reverse_3d;
-    P->fwd    = helmert_forward;
-    P->inv    = helmert_reverse;
 
     Q = (struct pj_opaque_helmert *)P->opaque;
 
diff --git a/test/gie/4D-API_cs2cs-style.gie b/test/gie/4D-API_cs2cs-style.gie
index bcdc256f..424501cc 100644
--- a/test/gie/4D-API_cs2cs-style.gie
+++ b/test/gie/4D-API_cs2cs-style.gie
@@ -185,7 +185,7 @@ operation proj=pipeline
                       k_0=1 x_0=2600000 y_0=1200000 ellps=bessel units=m
 -------------------------------------------------------------------------------
 tolerance 20 cm
-accept    7.438632495 46.951082877
+accept    7.438632495 46.951082877 0
 expect    2600000.0  1200000.0
 -------------------------------------------------------------------------------
 
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index 15bbaf05..32aad873 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -63,6 +63,11 @@ expect  -0.001790494 0.000895246
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+operation +proj=aea +ellps=GRS80 +lat_1=900
+expect  failure errno lat_larger_than_90
+
+operation +proj=aea +ellps=GRS80 +lat_2=900
+expect  failure errno lat_larger_than_90
 
 ===============================================================================
 Azimuthal Equidistant
@@ -1850,7 +1855,7 @@ Hammer & Eckert-Greifendorff
 ===============================================================================
 
 -------------------------------------------------------------------------------
-operation +proj=hammer   +a=6400000    +lat_1=0.5 +lat_2=2
+operation +proj=hammer   +a=6400000
 -------------------------------------------------------------------------------
 tolerance 0.1 mm
 accept  2 1
@@ -1873,6 +1878,12 @@ accept  -200 -100
 expect  -0.001790493 -0.000895247
 
 
+-------------------------------------------------------------------------------
+operation +proj=hammer   +a=6400000 +W=1
+-------------------------------------------------------------------------------
+accept  -180 0
+expect  failure errno tolerance_condition
+
 ===============================================================================
 Hatano Asymmetrical Equal Area
 	PCyl, Sph.
@@ -4698,7 +4709,7 @@ Oblique Stereographic Alternative
 ===============================================================================
 
 -------------------------------------------------------------------------------
-operation +proj=sterea   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +n=0.5
+operation +proj=sterea   +ellps=GRS80
 -------------------------------------------------------------------------------
 tolerance 0.1 mm
 accept  2 1
@@ -4721,7 +4732,7 @@ accept  -200 -100
 expect  -0.001796631 -0.000904369
 
 -------------------------------------------------------------------------------
-operation +proj=sterea   +R=6400000    +lat_1=0.5 +lat_2=2 +n=0.5
+operation +proj=sterea   +R=6400000
 -------------------------------------------------------------------------------
 tolerance 0.1 mm
 accept  2 1
@@ -4743,6 +4754,58 @@ expect  -0.001790493 0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+-------------------------------------------------------------------------------
+operation +proj=sterea   +ellps=GRS80    +lat_0=90
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  0 90
+expect  0 0
+accept  0 89
+expect  0.000000000000   -111696.700323081997
+accept  0 45
+expect  0.000000000000   -5291160.727484324016
+accept  0 0
+expect  0.000000000000   -12713600.098641794175
+
+-------------------------------------------------------------------------------
+operation +proj=sterea   +ellps=GRS80    +lat_0=89
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  0 90
+expect  0.000000000000   111696.700314355621
+accept  0 89
+expect  0.000000000000   0
+accept  0 45
+expect  0.000000000000   -5160845.342319893651
+accept  0 0
+expect  0.000000000000   -12493602.143489977345
+
+-------------------------------------------------------------------------------
+operation +proj=sterea   +ellps=GRS80    +lat_0=-90
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  0 -90
+expect  0 0
+accept  0 -89
+expect  0.000000000000   111696.700323081997
+accept  0 -45
+expect  0.000000000000   5291160.727484324016
+accept  0 0
+expect  0.000000000000   12713600.098641794175
+
+-------------------------------------------------------------------------------
+operation +proj=sterea   +ellps=GRS80   +lat_0=-89
+-------------------------------------------------------------------------------
+tolerance 0.1 mm
+accept  0 -90
+expect  0.000000000000   -111696.700314355621
+accept  0 -89
+expect  0.000000000000   0
+accept  0 -45
+expect  0.000000000000   5160845.342319893651
+accept  0 0
+expect  0.000000000000   12493602.143489977345
+
 
 ===============================================================================
 Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)
diff --git a/test/gie/ellipsoid.gie b/test/gie/ellipsoid.gie
index 2a6d02a5..2a7bdd60 100644
--- a/test/gie/ellipsoid.gie
+++ b/test/gie/ellipsoid.gie
@@ -67,6 +67,9 @@ expect    failure errno major_axis_not_given
 operation proj=merc +R=0
 expect    failure errno major_axis_not_given
 
+operation +proj=merc +R_a +a=2 +f=2
+expect  failure errno major_axis_not_given
+
 operation
 expect    failure
 operation cobra
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index d1c978a3..e16443d4 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -433,6 +433,11 @@ expect     failure   errno invalid_arg
 operation  proj=helmert transpose
 expect     failure   errno invalid_arg
 
+# Use of 2D Helmert interface with 3D Helmert setup
+operation  +proj=ob_tran +o_proj=helmert +o_lat_p=0
+direction  inverse
+accept     0 0
+expect     failure   errno 22
 
 -------------------------------------------------------------------------------
 Molodensky-Badekas from IOGP Guidance 7.2, Transformation from La Canoa to REGVEN