Merge pull request #1322 from rouault/ossfuzz_fixes

Various ossfuzz fixes

7 files changed, 43 insertions, 19 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;