IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

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:
6 edited
2 copied

Legend:

Unmodified
Added
Removed

  • trunk

  • trunk/Ohana

  • trunk/Ohana/src/opihi/cmd.data/Makefile

    r39638 r39926  
    124124$(SRC)/resize.$(ARCH).o         \
    125125$(SRC)/relocate.$(ARCH).o       \
     126$(SRC)/rndseed.$(ARCH).o                \
    126127$(SRC)/roll.$(ARCH).o           \
    127128$(SRC)/rotate.$(ARCH).o \

  • trunk/Ohana/src/opihi/cmd.data/init.c

    r39638 r39926  
    113113int reindex          PROTO((int, char **));
    114114int relocate         PROTO((int, char **));
     115int rndseed          PROTO((int, char **));
    115116int roll             PROTO((int, char **));
    116117int rotate           PROTO((int, char **));
    … …  
    293294  {1, "relocate",     relocate,         "set graphics/image window position"},
    294295  {1, "roll",         roll,             "roll image to new start point"},
     296  {1, "rndseed",      rndseed,          "set the pseudo-random seed"},
    295297  {1, "rotate",       rotate,           "rotate image"},
    296298  {1, "save",         save,             "save an SAOimage style image overlay"},

  • trunk/Ohana/src/opihi/cmd.data/match2d.c

    r36084 r39926  
    11# include "data.h"
    22
    3 int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2);
     3int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch);
    44int find_matches2d_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index);
     5
     6int find_matches2d_sphere (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch);
     7int find_matches2d_sphere_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index);
    58
    69// match2d (X1) (Y1) (X2) (Y2) (Radius) [-index1 (index1)] [-index2 (index2)] [-nomatch1 nomatch1] [-nomatch2 nomatch2]
    … …  
    2326  }
    2427
     28  int SPHERE_DISTANCE = FALSE;
     29  if ((N = get_argument (argc, argv, "-sphere"))) {
     30    remove_argument (N, &argc, argv);
     31    SPHERE_DISTANCE = TRUE;
     32  }
     33  if ((N = get_argument (argc, argv, "-sky"))) {
     34    remove_argument (N, &argc, argv);
     35    SPHERE_DISTANCE = TRUE;
     36  }
     37
    2538  if ((N = get_argument (argc, argv, "-index1"))) {
    2639    remove_argument (N, &argc, argv);
    … …  
    3750  } else {
    3851    if ((index2 = SelectVector ("index2", ANYVECTOR, TRUE)) == NULL) return (FALSE);   
     52  }
     53
     54 
     55  Vector *radiusMatch = NULL;
     56  if ((N = get_argument (argc, argv, "-radius"))) {
     57    if (CLOSEST) {
     58      gprint (GP_ERR, "error: -radius and -closest are currently incompatible\n");
     59      return (FALSE);
     60    }
     61    remove_argument (N, &argc, argv);
     62    if ((radiusMatch = SelectVector (argv[N], ANYVECTOR, TRUE)) == NULL) return (FALSE);   
     63    remove_argument (N, &argc, argv);
    3964  }
    4065
    … …  
    84109  }
    85110
    86   if (CLOSEST) {
     111  if (SPHERE_DISTANCE) {
     112    if (CLOSEST) {
     113      find_matches2d_sphere_closest (X1vec, Y1vec, X2vec, Y2vec, Radius, index1);
     114      find_matches2d_sphere_closest (X2vec, Y2vec, X1vec, Y1vec, Radius, index2);
     115    } else {
     116      find_matches2d_sphere (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2, radiusMatch);
     117    }
     118  } else {
     119    if (CLOSEST) {
    87120      find_matches2d_closest (X1vec, Y1vec, X2vec, Y2vec, Radius, index1);
    88121      find_matches2d_closest (X2vec, Y2vec, X1vec, Y1vec, Radius, index2);
    89   } else {
    90       find_matches2d (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2);
    91   }
    92 
     122    } else {
     123      find_matches2d (X1vec, Y1vec, X2vec, Y2vec, Radius, index1, index2, radiusMatch);
     124    }
     125  }
    93126  return (TRUE);
    94127
    … …  
    108141  gprint (GP_ERR, "use 'reindex' to generate new vectors based on these index vectors\n");
    109142
     143  gprint (GP_ERR, "if -sphere or -sky is supplied, (x1,y1) and (x2,y2) are treaded as (ra,dec) or (long,lat) pairs in degrees\n");
     144
     145  gprint (GP_ERR, "if -radius (vector) is supplied, the vector will be filled with the distance between the matched pairs\n");
     146  gprint (GP_ERR, "  not valid with -closest\n");
    110147  return FALSE;
    111148}
    112149
    113150// we are not defining a relative offset DX,DY for now
    114 int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2) {
     151int find_matches2d (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch) {
    115152 
    116153  off_t i, j, first_j, I, J, *N1, *N2, Nmatch, NMATCH, DMATCH;
    … …  
    122159  ResetVector (index1, OPIHI_INT, NMATCH);
    123160  ResetVector (index2, OPIHI_INT, NMATCH);
     161  if (radiusMatch) ResetVector (radiusMatch, OPIHI_FLT, NMATCH);
    124162
    125163  ALLOCATE (N1, off_t, X1->Nelements);
    … …  
    154192        index1->elements.Int[Nmatch] = I;
    155193        index2->elements.Int[Nmatch] = J;
     194        if (radiusMatch) radiusMatch->elements.Flt[Nmatch] = dR;
    156195
    157196        // XXX track matches 1 and 2 with internal vector, save new nomatch index vectors
    … …  
    163202          REALLOCATE (index1->elements.Int, opihi_int, NMATCH);
    164203          REALLOCATE (index2->elements.Int, opihi_int, NMATCH);
     204          if (radiusMatch) { REALLOCATE (radiusMatch->elements.Flt, opihi_flt, NMATCH); }
    165205        }
    166206      }
    … …  
    171211  index1->Nelements = Nmatch;
    172212  index2->Nelements = Nmatch;
     213  if (radiusMatch) radiusMatch->Nelements = Nmatch;
    173214
    174215  free (N1);
    … …  
    246287  return (TRUE);
    247288}
     289
     290double gcdist (double r1, double d1, double r2, double d2) {
     291  double num,den;
     292  r1 *= RAD_DEG;
     293  d1 *= RAD_DEG;
     294  r2 *= RAD_DEG;
     295  d2 *= RAD_DEG;
     296
     297  num = sqrt(pow((cos(d2) * sin(r2 - r1)),2) +
     298             pow((cos(d1) * sin(d2) -
     299                  sin(d1) * cos(d2) * cos(r2 - r1)),2));
     300  den = (sin(d1) * sin(d2) + cos(d1) * cos(d2) * cos(r2 - r1));
     301  return(atan2(num,den) * (180 / M_PI));
     302}
     303
     304typedef struct {
     305  double sD;
     306  double cD;
     307  double sR;
     308  double cR;
     309} Match2D_PreCalc;
     310
     311double gcdist_PreCalc_v0 (Match2D_PreCalc *p1, Match2D_PreCalc *p2) {
     312  double num,den;
     313
     314  num = sqrt(pow((p2->cD * (p2->sR*p1->cR - p1->sR*p2->cR)),2) +
     315             pow((p1->cD * p2->sD -
     316                  p1->sD * p2->cD * (p2->cR*p1->cR + p2->sR*p1->sR)),2));
     317  den = (p1->sD * p2->sD + p1->cD * p2->cD * (p2->cR*p1->cR + p2->sR*p1->sR));
     318  return(atan2(num,den) * (180 / M_PI));
     319}
     320
     321// we are not defining a relative offset DX,DY for now
     322double gcdist_PreCalc (Match2D_PreCalc *p1, Match2D_PreCalc *p2) {
     323  double num,den;
     324
     325  // double Qd = p1->sD * p1->cD * p2->sD * p2->cD;
     326  // double Qr = p1->sR * p1->cR * p2->sR * p2->cR;
     327
     328  //  double Xa
     329  //    = SQ(p2->cD * p1->cR * p2->sR)
     330  //    + SQ(p2->cD * p1->sR * p2->cR) -
     331  //    - 2 * SQ(p2->cD) * Qr;
     332
     333  double cdR = (p2->cR*p1->cR + p2->sR*p1->sR);
     334
     335  double Xa = SQ(p2->cD * p2->sR * p1->cR - p2->cD * p1->sR * p2->cR);
     336  double Xb = SQ(p1->cD * p2->sD - p1->sD * p2->cD * cdR);
     337
     338  num = sqrt(Xa + Xb);
     339  den = (p1->sD * p2->sD + p1->cD * p2->cD * cdR);
     340  return(atan2(num,den) * (180 / M_PI));
     341}
     342
     343// we are not defining a relative offset DX,DY for now
     344int find_matches2d_sphere (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index1, Vector *index2, Vector *radiusMatch) {
     345 
     346  off_t i, j, first_j, I, J, *N1, *N2, Nmatch, NMATCH, DMATCH;
     347  double dY, dR;
     348
     349  NMATCH = MAX(MAX(0.05*X1->Nelements, 0.05*X2->Nelements), 1000);
     350  DMATCH = NMATCH;
     351
     352  ResetVector (index1, OPIHI_INT, NMATCH);
     353  ResetVector (index2, OPIHI_INT, NMATCH);
     354  if (radiusMatch) ResetVector (radiusMatch, OPIHI_FLT, NMATCH);
     355
     356  ALLOCATE (N1, off_t, X1->Nelements);
     357  ALLOCATE (N2, off_t, X2->Nelements);
     358
     359  ALLOCATE_PTR (A1, Match2D_PreCalc, X1->Nelements);
     360  ALLOCATE_PTR (A2, Match2D_PreCalc, X2->Nelements);
     361
     362  for (i = 0; i < X1->Nelements; i++) {
     363    A1[i].sR = sin(RAD_DEG*X1->elements.Flt[i]);
     364    A1[i].cR = cos(RAD_DEG*X1->elements.Flt[i]);
     365    A1[i].sD = sin(RAD_DEG*Y1->elements.Flt[i]);
     366    A1[i].cD = cos(RAD_DEG*Y1->elements.Flt[i]);
     367    N1[i] = i;
     368  }
     369  for (i = 0; i < X2->Nelements; i++) {
     370    A2[i].sR = sin(RAD_DEG*X2->elements.Flt[i]);
     371    A2[i].cR = cos(RAD_DEG*X2->elements.Flt[i]);
     372    A2[i].sD = sin(RAD_DEG*Y2->elements.Flt[i]);
     373    A2[i].cD = cos(RAD_DEG*Y2->elements.Flt[i]);
     374    N2[i] = i;
     375  }
     376
     377  // sort from one pole to the other
     378  sort_coords_indexonly (Y1->elements.Flt, X1->elements.Flt, N1, X1->Nelements);
     379  sort_coords_indexonly (Y2->elements.Flt, X2->elements.Flt, N2, X2->Nelements);
     380
     381  Nmatch = 0;
     382  for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) {
     383    I = N1[i];
     384    J = N2[j];
     385
     386    // we can use dY as minimal requirement: if dY > Radius, we are too far apart
     387    dY = Y1->elements.Flt[I] - Y2->elements.Flt[J];
     388
     389    if (dY <= -1.02*Radius) { i++; continue; }
     390    if (dY >= +1.02*Radius) { j++; continue; }
     391
     392    // look for all matches of list2() to list1(i)
     393    first_j = j;
     394    for (j = first_j; (dY > -1.02*Radius) && (j < X2->Nelements); j++) {
     395      J = N2[j];
     396
     397      dR = gcdist_PreCalc (&A1[I], &A2[J]);
     398      // dR = gcdist (X1->elements.Flt[I], Y1->elements.Flt[I], X2->elements.Flt[J], Y2->elements.Flt[J]);
     399
     400      if (dR < Radius) {
     401        index1->elements.Int[Nmatch] = I;
     402        index2->elements.Int[Nmatch] = J;
     403        if (radiusMatch) radiusMatch->elements.Flt[Nmatch] = dR;
     404
     405        // XXX track matches 1 and 2 with internal vector, save new nomatch index vectors
     406        // after this loop
     407
     408        Nmatch ++;
     409        if (Nmatch >= NMATCH) {
     410          NMATCH += DMATCH;
     411          REALLOCATE (index1->elements.Int, opihi_int, NMATCH);
     412          REALLOCATE (index2->elements.Int, opihi_int, NMATCH);
     413          if (radiusMatch) { REALLOCATE (radiusMatch->elements.Flt, opihi_flt, NMATCH); }
     414        }
     415      }
     416    }
     417    j = first_j;
     418    i++;
     419  }
     420  index1->Nelements = Nmatch;
     421  index2->Nelements = Nmatch;
     422  if (radiusMatch) radiusMatch->Nelements = Nmatch;
     423
     424  free (A1);
     425  free (A2);
     426
     427  free (N1);
     428  free (N2);
     429
     430  return (TRUE);
     431}
     432
     433// we are not defining a relative offset DX,DY for now
     434int find_matches2d_sphere_closest (Vector *X1, Vector *Y1, Vector *X2, Vector *Y2, double Radius, Vector *index) {
     435 
     436  off_t i, j, Jmin, Ji, I, J, *N1, *N2;
     437  double dY, dR, Rmin;
     438
     439  ResetVector (index, OPIHI_INT, X1->Nelements);
     440  for (i = 0; i < index->Nelements; i++) { index->elements.Int[i] = -1; }
     441
     442  ALLOCATE (N1, off_t, X1->Nelements);
     443  ALLOCATE (N2, off_t, X2->Nelements);
     444
     445  ALLOCATE_PTR (A1, Match2D_PreCalc, X1->Nelements);
     446  ALLOCATE_PTR (A2, Match2D_PreCalc, X2->Nelements);
     447
     448  for (i = 0; i < X1->Nelements; i++) {
     449    A1[i].sR = sin(RAD_DEG*X1->elements.Flt[i]);
     450    A1[i].cR = cos(RAD_DEG*X1->elements.Flt[i]);
     451    A1[i].sD = sin(RAD_DEG*Y1->elements.Flt[i]);
     452    A1[i].cD = cos(RAD_DEG*Y1->elements.Flt[i]);
     453    N1[i] = i;
     454  }
     455  for (i = 0; i < X2->Nelements; i++) {
     456    A2[i].sR = sin(RAD_DEG*X2->elements.Flt[i]);
     457    A2[i].cR = cos(RAD_DEG*X2->elements.Flt[i]);
     458    A2[i].sD = sin(RAD_DEG*Y2->elements.Flt[i]);
     459    A2[i].cD = cos(RAD_DEG*Y2->elements.Flt[i]);
     460    N2[i] = i;
     461  }
     462
     463  // sort from one pole to the other
     464  sort_coords_indexonly (Y1->elements.Flt, X1->elements.Flt, N1, X1->Nelements);
     465  sort_coords_indexonly (Y2->elements.Flt, X2->elements.Flt, N2, X2->Nelements);
     466
     467  for (i = j = 0; (i < X1->Nelements) && (j < X2->Nelements);) {
     468    I = N1[i];
     469    J = N2[j];
     470
     471    // we can use dY as minimal requirement: if dY > Radius, we are too far apart
     472    dY = Y1->elements.Flt[I] - Y2->elements.Flt[J];
     473
     474    if (dY <= -1.02*Radius) {
     475      // no match in list 2 to this entry
     476      index->elements.Int[I] = -1; // (probably not needed --- index is init'ed above)o
     477      i++;
     478      continue;
     479    }
     480    if (dY >= +1.02*Radius) { j++; continue; }
     481
     482    // look for all matches of list2() to list1(i)
     483    Jmin = -1;
     484    Rmin = Radius;
     485    for (Ji = j; (dY > -1.02*Radius) && (Ji < X2->Nelements); Ji++) {
     486      J = N2[Ji];
     487
     488      dR = gcdist_PreCalc (&A1[I], &A2[J]);
     489      // dR = gcdist (X1->elements.Flt[I], Y1->elements.Flt[I], X2->elements.Flt[J], Y2->elements.Flt[J]);
     490
     491      if (dR < Rmin) {
     492        Rmin = dR;
     493        Jmin = J;
     494      }
     495    }
     496   
     497
     498    // no match in list 2 to this entry
     499    if (Jmin == -1) {
     500      index->elements.Int[I] = -1;
     501      i++;
     502      continue;
     503    }
     504    index->elements.Int[I] = Jmin;
     505    i++;
     506  }
     507
     508  free (A1);
     509  free (A2);
     510
     511  free (N1);
     512  free (N2);
     513
     514  return (TRUE);
     515}
     516
     517

  • trunk/Ohana/src/opihi/cmd.data/stats-new.c

    r15274 r39926  
    129129  ALLOCATE (values, float, Nsample);
    130130   
    131   A = time(NULL);
    132   for (B = 0; A == time(NULL); B++);
    133   srand48(B);
     131  // srand48() is called by startup.c
    134132 
    135133  *buffer = (float *) matrix[0].buffer;
Note: See TracChangeset for help on using the changeset viewer.