IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 39926 for trunk/psLib


Ignore:
Timestamp:
Jan 6, 2017, 11:30:10 AM (10 years ago)
Author:
eugene
Message:

merging changes from czw dev branch (compare with r39924)

Location:
trunk
Files:
5 edited

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/psLib/src/astro/psCoord.c

    r19504 r39926  
    234234    const psPlane* coords)
    235235{
     236  if (!transform) {
     237    fprintf (stderr, "problem\n");
     238  }
    236239    PS_ASSERT_PTR_NON_NULL(transform, NULL);
    237240    PS_ASSERT_PTR_NON_NULL(transform->x, NULL);
    … …  
    305308void projectionFree(psProjection *p)
    306309{
    307     // There are no dynamically allocated items
     310  if (!p) return;
     311  if (!p->radial) return;
     312  psFree (p->radial);
    308313}
    309314
    … …  
    321326    p->Ys = Ys;
    322327    p->type = type;
     328    p->radial = NULL;
    323329
    324330    psMemSetDeallocator(p, (psFreeFunc) projectionFree);
    … …  
    368374      case PS_PROJ_ZEA:
    369375      case PS_PROJ_ZPL:
     376      case PS_PROJ_ZPN:
    370377        class = PS_PROJECTION_CLASS_ZENITHAL;
    371378        break;
    … …  
    442449            out->y = -zeta * cosPhiCT;
    443450            break;
     451
     452          case PS_PROJ_ZPN: {
     453            // the forward projection is:
     454            // theta = atan2(stht, ctht)
     455            // gamma = (pi/2 - theta) : theta in radians
     456            // Ro = sum (P_i gamma^i)
     457            // R  = (180/pi) Ro
     458           
     459            // Ro = (pi/180)(90 - theta)
     460            // R = (180/pi)sum (P_i R^i)
     461           
     462            // is ZPN defined for Npolyterms = 0 or 1?
     463
     464            double cosTheta = hypot(sinPhiCT, cosPhiCT);
     465            double theta = atan2 (sinTheta, cosTheta);
     466           
     467            double gamma = M_PI_2 - theta;
     468           
     469            // i = 0 .. N - 1 (1 <= Npolyterms <= 21)
     470            double Ro = 0.0;
     471            for (int i = projection->radial->n - 1; i > 0; i--) {
     472              double Pi = projection->radial->data.F64[i];
     473              Ro = (Ro + Pi)*gamma;
     474            }
     475            Ro += projection->radial->data.F64[0];
     476           
     477            out->x = (cosTheta == 0.0) ? 0.0 : +Ro * sinPhiCT / cosTheta ;
     478            out->y = (cosTheta == 0.0) ? 0.0 : -Ro * cosPhiCT / cosTheta ;
     479            break;
     480          }
     481
    444482          default:
    445483            psAbort("invalid projection");
    … …  
    525563      case PS_PROJ_ZEA:
    526564      case PS_PROJ_ZPL:
     565      case PS_PROJ_ZPN:
    527566        class = PS_PROJECTION_CLASS_ZENITHAL;
    528567        break;
    … …  
    574613          case PS_PROJ_SIN:
    575614            // Orhtographic deprojection
     615            cosTheta = R;
    576616            sinTheta = sqrt (1 - R*R);
    577             cosTheta = R;
    578617            break;
    579618          case PS_PROJ_STG:
    580             psAbort("STG not defined");
     619            sinTheta = (4 - R) / (4 + R);
     620            cosTheta = sqrt (1 - sinTheta*sinTheta);
    581621            break;
    582622          case PS_PROJ_ZEA:
    … …  
    587627            cosTheta = sqrt (1 - PS_SQR(sinTheta));
    588628            break;
     629
     630      case PS_PROJ_ZPN:
     631
     632        // the forward projection is:
     633        // theta = atan2(stht, ctht)
     634        // gamma = (pi/2 - theta) : theta in radians
     635        // Ro = sum (P_i gamma^i)
     636        // R  = (180/pi) Ro
     637
     638        // given R, we need to find theta:
     639        // Ro = R * (pi / 180) = sum (P_i gamma^i)
     640        // solve sum (P_i gamma^i) - Ro = 0 using Newton-Raphson
     641
     642        // use Ro to get a guess for gamma and iterate
     643
     644        {
     645          // find the roots of f(gamma) - R = 0
     646          // starting guess for gamma is (R - P0) / P1
     647          double gamma = (R - projection->radial->data.F64[0]) / projection->radial->data.F64[1];
     648         
     649          for (int iter = 0; iter < 5; iter++) {
     650           
     651            double Rc = 0.0; // this will hold the ander
     652            double dR = 0.0;
     653            for (int i = projection->radial->n - 1; i > 1; i--) {
     654              double Pi = projection->radial->data.F64[i];
     655              Rc = (Rc + Pi)*gamma;
     656              dR = (dR + i*Pi)*gamma;
     657            }
     658            double P0 = projection->radial->data.F64[0];
     659            double P1 = projection->radial->data.F64[1];
     660            Rc = (Rc + P1)*gamma + P0;
     661            dR = (dR + P1);
     662
     663            double gamma_new = gamma - (Rc - R) / dR;
     664            gamma = gamma_new;
     665          }
     666         
     667          double theta = M_PI_2 - gamma ;
     668          cosTheta = cos (theta);
     669          sinTheta = sin (theta);
     670          break;
     671        }
     672       
    589673          default:
    590674            psAbort("invalid projection");
    … …  
    9491033        myPT = psPlaneTransformAlloc(order, order);
    9501034    } else {
     1035      // the user has supplied a model with a specific order : fit that order
     1036      myPT = psMemIncrRefCounter(out); // we need to return something which can be freed
     1037# if (0)
    9511038        if ((out->x->nX == order) && (out->x->nY == order) &&
    9521039            (out->y->nX == order) && (out->y->nX == order)) {
    … …  
    9561043            myPT = psPlaneTransformAlloc(order, order);
    9571044        }
     1045# endif
    9581046    }
    9591047
    … …  
    11831271
    11841272    switch (type) {
    1185       case PS_PROJ_LIN:
    1186         psStringAppend (&name, "%s-LIN", prefix);
    1187         return name;
    1188       case PS_PROJ_PLY:
    1189         psStringAppend (&name, "%s-PLY", prefix);
    1190         return name;
    1191       case PS_PROJ_WRP:
    1192         psStringAppend (&name, "%s-WRP", prefix);
    1193         return name;
    1194       case PS_PROJ_TAN:
    1195         psStringAppend (&name, "%s-TAN", prefix);
    1196         return name;
    1197       case PS_PROJ_TNX:
    1198         psStringAppend (&name, "%s-TNX", prefix);
    1199         return name;
    1200       case PS_PROJ_DIS:
    1201         psStringAppend (&name, "%s-DIS", prefix);
    1202         return name;
    1203       case PS_PROJ_SIN:
    1204         psStringAppend (&name, "%s-SIN", prefix);
    1205         return name;
    1206       case PS_PROJ_STG:
    1207         psStringAppend (&name, "%s-STG", prefix);
    1208         return name;
    1209       case PS_PROJ_AIT:
    1210         psStringAppend (&name, "%s-AIT", prefix);
    1211         return name;
    1212       case PS_PROJ_PAR:
    1213         psStringAppend (&name, "%s-PAR", prefix);
    1214         return name;
    1215       case PS_PROJ_GLS:
    1216         psStringAppend (&name, "%s-GLS", prefix);
    1217         return name;
    1218       case PS_PROJ_CAR:
    1219         psStringAppend (&name, "%s-CAR", prefix);
    1220         return name;
    1221       case PS_PROJ_MER:
    1222         psStringAppend (&name, "%s-MER", prefix);
    1223         return name;
     1273      case PS_PROJ_LIN: psStringAppend (&name, "%s-LIN", prefix); return name;
     1274      case PS_PROJ_PLY: psStringAppend (&name, "%s-PLY", prefix); return name;
     1275      case PS_PROJ_WRP: psStringAppend (&name, "%s-WRP", prefix); return name;
     1276      case PS_PROJ_TAN: psStringAppend (&name, "%s-TAN", prefix); return name;
     1277
     1278      case PS_PROJ_DIS: psStringAppend (&name, "%s-DIS", prefix); return name;
     1279      case PS_PROJ_SIN: psStringAppend (&name, "%s-SIN", prefix); return name;
     1280      case PS_PROJ_STG: psStringAppend (&name, "%s-STG", prefix); return name;
     1281      case PS_PROJ_TNX: psStringAppend (&name, "%s-TNX", prefix); return name;
     1282
     1283      case PS_PROJ_ZEA: psStringAppend (&name, "%s-ZEA", prefix); return name;
     1284      case PS_PROJ_ZPL: psStringAppend (&name, "%s-ZPL", prefix); return name;
     1285      case PS_PROJ_ZPN: psStringAppend (&name, "%s-ZPN", prefix); return name;
     1286      case PS_PROJ_AIT: psStringAppend (&name, "%s-AIT", prefix); return name;
     1287
     1288      case PS_PROJ_PAR: psStringAppend (&name, "%s-PAR", prefix); return name;
     1289      case PS_PROJ_GLS: psStringAppend (&name, "%s-GLS", prefix); return name;
     1290      case PS_PROJ_CAR: psStringAppend (&name, "%s-CAR", prefix); return name;
     1291      case PS_PROJ_MER: psStringAppend (&name, "%s-MER", prefix); return name;
    12241292
    12251293      default:
    … …  
    12381306    if (!strcmp (&name[4], "-WRP")) return PS_PROJ_WRP;
    12391307    if (!strcmp (&name[4], "-TAN")) return PS_PROJ_TAN;
    1240     if (!strcmp (&name[4], "-TNX")) return PS_PROJ_TNX;
    12411308    if (!strcmp (&name[4], "-DIS")) return PS_PROJ_DIS;
    12421309    if (!strcmp (&name[4], "-SIN")) return PS_PROJ_SIN;
    12431310    if (!strcmp (&name[4], "-STG")) return PS_PROJ_STG;
     1311    if (!strcmp (&name[4], "-TNX")) return PS_PROJ_TNX;
     1312    if (!strcmp (&name[4], "-ZEA")) return PS_PROJ_ZEA;
     1313    if (!strcmp (&name[4], "-ZPL")) return PS_PROJ_ZPL;
     1314    if (!strcmp (&name[4], "-ZPN")) return PS_PROJ_ZPN;
    12441315    if (!strcmp (&name[4], "-AIT")) return PS_PROJ_AIT;
    12451316    if (!strcmp (&name[4], "-PAR")) return PS_PROJ_PAR;

  • trunk/psLib/src/astro/psCoord.h

    r14452 r39926  
    130130    PS_PROJ_ZEA,                        ///< Sine projection
    131131    PS_PROJ_ZPL,                        ///< Sine projection
     132    PS_PROJ_ZPN,                        ///< Sine projection
    132133    PS_PROJ_AIT,                        ///< Aitoff projection
    133134    PS_PROJ_PAR,                        ///< Par projection
    … …  
    150151    double Ys;                         ///< plate-scale in Y direction
    151152    psProjectionType type;             ///< Projection type
    152 }
    153 psProjection;
     153    psVector *radial;          ///< radial distortion terms
     154} psProjection;
    154155
    155156/** Allocates a psPlane

  • trunk/psLib/src/fits/psFitsScale.c

    r39509 r39926  
    189189                // Desperate retry
    190190                mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT);
    191                 stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
     191                if ((STDEV_STAT == PS_STAT_SAMPLE_QUARTILE)||(STDEV_STAT == PS_STAT_ROBUST_QUARTILE)) {
     192                  stdev = 0.7415 * (stats->sampleUQ - stats->sampleLQ);
     193                }
     194                else {
     195                  stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
     196                }
    192197            }
    193198        } else {
    194199            // Retry with all available pixels
    195200            mean = psStatsGetValue(stats, MEAN_STAT);
    196             stdev = psStatsGetValue(stats, STDEV_STAT);
     201            if ((STDEV_STAT == PS_STAT_SAMPLE_QUARTILE)||(STDEV_STAT == PS_STAT_ROBUST_QUARTILE)) {
     202              stdev = 0.7415 * (stats->sampleUQ - stats->sampleLQ);
     203            }
     204            else {
     205              stdev = psStatsGetValue(stats, STDEV_STAT);
     206            }
    197207        }
    198208    } else {
    199209        // First attempt
    200210        mean = psStatsGetValue(stats, MEAN_STAT);
    201         stdev = psStatsGetValue(stats, STDEV_STAT);
     211        if ((STDEV_STAT == PS_STAT_SAMPLE_QUARTILE)||(STDEV_STAT == PS_STAT_ROBUST_QUARTILE)) {
     212          stdev = 0.7415 * (stats->sampleUQ - stats->sampleLQ);
     213        }
     214        else {
     215          stdev = psStatsGetValue(stats, STDEV_STAT);
     216        }
    202217    }
    203218    psFree(rng);
    … …  
    554569                // Desperate retry
    555570                mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT);
    556                 stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
    557             }
     571                if ((DESPERATE_STDEV_STAT == PS_STAT_SAMPLE_QUARTILE)||(DESPERATE_STDEV_STAT == PS_STAT_ROBUST_QUARTILE)) {
     572                  stdev = 0.7415 * (stats->sampleUQ - stats->sampleLQ); // Calculate a "sigma" based on the interquartile distance.
     573                }
     574                else {
     575                 
     576                  stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
     577                }
     578            }   
    558579        } else {
    559580            // Retry with all available pixels
    560581            mean = psStatsGetValue(stats, MEAN_STAT);
    561             stdev = psStatsGetValue(stats, STDEV_STAT);
     582            if ((STDEV_STAT == PS_STAT_SAMPLE_QUARTILE)||(STDEV_STAT == PS_STAT_ROBUST_QUARTILE)) {
     583              stdev = 0.7415 * (stats->sampleUQ - stats->sampleLQ); // Calculate a "sigma" based on the interquartile distance.
     584            }
     585            else {
     586              stdev = psStatsGetValue(stats, STDEV_STAT);
     587            }
    562588        }
    563589    } else {
    564590        // First attempt
    565591        mean = psStatsGetValue(stats, MEAN_STAT);
    566         stdev = psStatsGetValue(stats, STDEV_STAT);
     592        if ((STDEV_STAT == PS_STAT_SAMPLE_QUARTILE)||(STDEV_STAT == PS_STAT_ROBUST_QUARTILE)) {
     593          stdev = 0.7415 * (stats->sampleUQ - stats->sampleLQ); // Calculate a "sigma" based on the interquartile distance.
     594        }
     595        else {
     596          stdev = psStatsGetValue(stats, STDEV_STAT);
     597        }
    567598    }
    568599    psFree(rng);
    … …  
    723754                // Desperate retry
    724755                mean = psStatsGetValue(stats, DESPERATE_MEAN_STAT);
    725                 stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
     756                if ((DESPERATE_STDEV_STAT == PS_STAT_SAMPLE_QUARTILE)||(DESPERATE_STDEV_STAT == PS_STAT_ROBUST_QUARTILE)) {
     757                  stdev = 0.7415 * (stats->sampleUQ - stats->sampleLQ); // Calculate a "sigma" based on the interquartile distance.
     758                }
     759                else {
     760                  stdev = psStatsGetValue(stats, DESPERATE_STDEV_STAT);
     761                }
    726762            }
    727763        } else {
    728764            // Retry with all available pixels
    729765            mean = psStatsGetValue(stats, MEAN_STAT);
    730             stdev = psStatsGetValue(stats, STDEV_STAT);
     766            if ((STDEV_STAT == PS_STAT_SAMPLE_QUARTILE)||(STDEV_STAT == PS_STAT_ROBUST_QUARTILE)) {
     767              stdev = 0.7415 * (stats->sampleUQ - stats->sampleLQ); // Calculate a "sigma" based on the interquartile distance.
     768            }
     769            else {
     770              stdev = psStatsGetValue(stats, STDEV_STAT);
     771            }
    731772        }
    732773    } else {
    733774        // First attempt
    734775        mean = psStatsGetValue(stats, MEAN_STAT);
    735         stdev = psStatsGetValue(stats, STDEV_STAT);
     776        if ((STDEV_STAT == PS_STAT_SAMPLE_QUARTILE)||(STDEV_STAT == PS_STAT_ROBUST_QUARTILE)) {
     777          stdev = 0.7415 * (stats->sampleUQ - stats->sampleLQ); // Calculate a "sigma" based on the interquartile distance.
     778        }
     779        else {
     780          stdev = psStatsGetValue(stats, STDEV_STAT);
     781        }
    736782    }
    737783    psFree(rng);

  • trunk/psLib/src/types/psMetadata.c

    r39598 r39926  
    13901390            break;
    13911391        case PS_DATA_F64:
    1392             fprintf(fd, "%f\n", item->data.F64);
     1392            fprintf(fd, "%g\n", item->data.F64);
    13931393            break;
    13941394        case PS_DATA_METADATA:
Note: See TracChangeset for help on using the changeset viewer.