← Back to branch summary

~ubuntu-branches/ubuntu/utopic/wcslib/utopic

« back to all changes in this revision

Viewing changes to C/wcs.c

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-05-13 10:29:41 UTC
  • mfrom: (15.1.3 sid)
  • Revision ID: package-import@ubuntu.com-20140513102941-c8643128v6rews10
Tags: 4.23-1
New upstream version

Show diffs side-by-side

added added

removed removed

Lines of Context:
C/wcs.c
1
1
/*============================================================================
2
2
 
3
 
  WCSLIB 4.22 - an implementation of the FITS WCS standard.
 
3
  WCSLIB 4.23 - an implementation of the FITS WCS standard.
4
4
  Copyright (C) 1995-2014, Mark Calabretta
5
5
 
6
6
  This file is part of WCSLIB.
22
22
 
23
23
  Author: Mark Calabretta, Australia Telescope National Facility, CSIRO.
24
24
  http://www.atnf.csiro.au/people/Mark.Calabretta
25
 
  $Id: wcs.c,v 4.22 2014/04/12 15:03:52 mcalabre Exp $
 
25
  $Id: wcs.c,v 4.23 2014/05/13 05:50:51 mcalabre Exp $
26
26
*===========================================================================*/
27
27
 
28
28
#include <math.h>
515
515
  memset(wcs->dateavg, 0, 72);
516
516
  wcs->mjdobs     = UNDEFINED;
517
517
  wcs->mjdavg     = UNDEFINED;
 
518
  wcs->velangl    = UNDEFINED;
518
519
 
519
520
  wcs->ntab = 0;
520
521
  wcs->tab  = 0x0;
548
549
  static const char *function = "wcssub";
549
550
 
550
551
  char *c, ctypei[16];
551
 
  int  axis, cubeface, dealloc, dummy, i, itab, j, k, latitude, longitude, m,
552
 
       *map = 0x0, msub, naxis, npv, nps, other, spectral, status, stokes;
 
552
  int  axis, cubeface, dealloc, dummy, i, itab, *itmp = 0x0, j, k, latitude,
 
553
       longitude, m, *map, msub, naxis, npv, nps, ntmp, other, spectral,
 
554
       status, stokes;
553
555
  const double *srcp;
554
556
  double *dstp;
555
557
  struct tabprm *tabp;
564
566
      "naxis must be positive (got %d)", naxis);
565
567
  }
566
568
 
567
 
  if (!(map = calloc(naxis, sizeof(int)))) {
568
 
    return wcserr_set(WCS_ERRMSG(WCSERR_MEMORY));
569
 
  }
570
 
 
571
569
  if (nsub == 0x0) {
572
570
    nsub = &dummy;
573
571
    *nsub = naxis;
575
573
    *nsub = naxis;
576
574
  }
577
575
 
 
576
  /* Allocate enough temporary storage to hold either axes[] or map[].*/
 
577
  ntmp = (*nsub <= naxis) ? naxis : *nsub;
 
578
  if (!(itmp = calloc(ntmp, sizeof(int)))) {
 
579
    return wcserr_set(WCS_ERRMSG(WCSERR_MEMORY));
 
580
  }
 
581
 
578
582
  if ((dealloc = (axes == 0x0))) {
579
583
    /* Construct an index array. */
580
584
    if (!(axes = calloc(naxis, sizeof(int)))) {
581
 
      free(map);
 
585
      free(itmp);
582
586
      return wcserr_set(WCS_ERRMSG(WCSERR_MEMORY));
583
587
    }
584
588
 
678
682
 
679
683
        /* This axis is wanted, but has it already been added? */
680
684
        for (k = 0; k < msub; k++) {
681
 
          if (map[k] == i+1) {
 
685
          if (itmp[k] == i+1) {
682
686
            break;
683
687
          }
684
688
        }
685
 
        if (k == msub) map[msub++] = i+1;
 
689
        if (k == msub) itmp[msub++] = i+1;
686
690
      }
687
691
 
688
692
    } else if (0 < axis && axis <= naxis) {
689
693
      /* Check that the requested axis has not already been added. */
690
694
      for (k = 0; k < msub; k++) {
691
 
        if (map[k] == axis) {
 
695
        if (itmp[k] == axis) {
692
696
          break;
693
697
        }
694
698
      }
695
 
      if (k == msub) map[msub++] = axis;
 
699
      if (k == msub) itmp[msub++] = axis;
696
700
 
697
701
    } else if (axis == 0) {
698
702
      /* Graft on a new axis. */
699
 
      map[msub++] = 0;
 
703
      itmp[msub++] = 0;
700
704
 
701
705
    } else {
702
706
      status = wcserr_set(WCS_ERRMSG(WCSERR_BAD_SUBIMAGE));
710
714
  }
711
715
 
712
716
  for (i = 0; i < *nsub; i++) {
713
 
    axes[i] = map[i];
 
717
    axes[i] = itmp[i];
714
718
  }
715
719
 
716
720
 
719
723
     axes[i] == 0 means to create a new axis,
720
724
      map[i] == j means that input axis i+1 goes to output axis j,
721
725
      map[i] == 0 means that input axis i+1 is not used. */
 
726
  map = itmp;
722
727
  for (i = 0; i < naxis; i++) {
723
728
    map[i] = 0;
724
729
  }
730
735
  }
731
736
 
732
737
  /* Check that the subimage coordinate system is separable. */
733
 
  if (*nsub < naxis) {
734
 
    srcp = wcssrc->pc;
735
 
    for (i = 0; i < naxis; i++) {
736
 
      for (j = 0; j < naxis; j++) {
737
 
        if (*(srcp++) == 0.0 || j == i) continue;
 
738
  srcp = wcssrc->pc;
 
739
  for (i = 0; i < naxis; i++) {
 
740
    for (j = 0; j < naxis; j++) {
 
741
      if (*(srcp++) == 0.0 || j == i) continue;
738
742
 
739
 
        if ((map[i] == 0) != (map[j] == 0)) {
740
 
          status = wcserr_set(WCS_ERRMSG(WCSERR_NON_SEPARABLE));
741
 
          goto cleanup;
742
 
        }
 
743
      if ((map[i] == 0) != (map[j] == 0)) {
 
744
        status = wcserr_set(WCS_ERRMSG(WCSERR_NON_SEPARABLE));
 
745
        goto cleanup;
743
746
      }
744
747
    }
745
748
  }
965
968
 
966
969
 
967
970
cleanup:
968
 
  if (map) free(map);
 
971
  if (itmp) free(itmp);
969
972
  if (dealloc) {
970
973
    free(axes);
971
974
  }
977
980
 
978
981
/*--------------------------------------------------------------------------*/
979
982
 
 
983
int wcscompare(
 
984
  int cmp,
 
985
  const struct wcsprm *wcs1,
 
986
  const struct wcsprm *wcs2,
 
987
  int *equal)
 
988
 
 
989
{
 
990
  int i, j, naxis, naxis2;
 
991
  double diff;
 
992
  int tab_equal;
 
993
  int status;
 
994
 
 
995
  if (wcs1  == 0x0) return WCSERR_NULL_POINTER;
 
996
  if (wcs2  == 0x0) return WCSERR_NULL_POINTER;
 
997
  if (equal == 0x0) return WCSERR_NULL_POINTER;
 
998
 
 
999
  *equal = 0;
 
1000
 
 
1001
  if (wcs1->naxis != wcs2->naxis) {
 
1002
    return 0;
 
1003
  }
 
1004
 
 
1005
  naxis = wcs1->naxis;
 
1006
  naxis2 = wcs1->naxis*wcs1->naxis;
 
1007
 
 
1008
  if (cmp & WCSCOMPARE_CRPIX) {
 
1009
    /* Don't compare crpix. */
 
1010
  } else if (cmp & WCSCOMPARE_TILING) {
 
1011
    for (i = 0; i < naxis; ++i) {
 
1012
      diff = wcs1->crpix[i] - wcs2->crpix[i];
 
1013
      if ((double)(int)(diff) != diff) {
 
1014
        return 0;
 
1015
      }
 
1016
    }
 
1017
  } else {
 
1018
    if (!wcsutil_Eq(naxis, wcs1->crpix, wcs2->crpix)) {
 
1019
      return 0;
 
1020
    }
 
1021
  }
 
1022
 
 
1023
  if (!wcsutil_Eq(naxis2, wcs1->pc, wcs2->pc) ||
 
1024
      !wcsutil_Eq(naxis, wcs1->cdelt, wcs2->cdelt) ||
 
1025
      !wcsutil_Eq(naxis, wcs1->crval, wcs2->crval) ||
 
1026
      !wcsutil_strEq(naxis, wcs1->cunit, wcs2->cunit) ||
 
1027
      !wcsutil_strEq(naxis, wcs1->ctype, wcs2->ctype) ||
 
1028
      wcs1->lonpole != wcs2->lonpole ||
 
1029
      wcs1->latpole != wcs2->latpole ||
 
1030
      wcs1->restfrq != wcs2->restfrq ||
 
1031
      wcs1->restwav != wcs2->restwav ||
 
1032
      wcs1->npv != wcs2->npv ||
 
1033
      wcs1->nps != wcs2->nps) {
 
1034
    return 0;
 
1035
  }
 
1036
 
 
1037
  /* Compare pv cards, which may not be in the same order */
 
1038
  for (i = 0; i < wcs1->npv; ++i) {
 
1039
    for (j = 0; j < wcs2->npv; ++j) {
 
1040
      if (wcs1->pv[i].i == wcs2->pv[j].i &&
 
1041
          wcs1->pv[i].m == wcs2->pv[j].m) {
 
1042
        if (wcs1->pv[i].value != wcs2->pv[j].value) {
 
1043
          return 0;
 
1044
        }
 
1045
        break;
 
1046
      }
 
1047
    }
 
1048
    /* We didn't find a match, so they are not equal */
 
1049
    if (j == wcs2->npv) {
 
1050
      return 0;
 
1051
    }
 
1052
  }
 
1053
 
 
1054
  /* Compare ps cards, which may not be in the same order */
 
1055
  for (i = 0; i < wcs1->nps; ++i) {
 
1056
    for (j = 0; j < wcs2->nps; ++j) {
 
1057
      if (wcs1->ps[i].i == wcs2->ps[j].i &&
 
1058
          wcs1->ps[i].m == wcs2->ps[j].m) {
 
1059
        if (strncmp(wcs1->ps[i].value, wcs2->ps[j].value, 72)) {
 
1060
          return 0;
 
1061
        }
 
1062
        break;
 
1063
      }
 
1064
    }
 
1065
    /* We didn't find a match, so they are not equal */
 
1066
    if (j == wcs2->nps) {
 
1067
      return 0;
 
1068
    }
 
1069
  }
 
1070
 
 
1071
  if (wcs1->flag != WCSSET || wcs2->flag != WCSSET) {
 
1072
    if (!wcsutil_Eq(naxis2, wcs1->cd, wcs2->cd) ||
 
1073
        !wcsutil_Eq(naxis, wcs1->crota, wcs2->crota) ||
 
1074
        wcs1->altlin != wcs2->altlin ||
 
1075
        wcs1->velref != wcs2->velref) {
 
1076
      return 0;
 
1077
    }
 
1078
  }
 
1079
 
 
1080
  if (!(cmp & WCSCOMPARE_ANCILLARY)) {
 
1081
    if (strncmp(wcs1->alt, wcs2->alt, 4) ||
 
1082
        wcs1->colnum != wcs2->colnum ||
 
1083
        !wcsutil_intEq(naxis, wcs1->colax, wcs2->colax) ||
 
1084
        !wcsutil_strEq(naxis, wcs1->cname, wcs2->cname) ||
 
1085
        !wcsutil_Eq(naxis, wcs1->crder, wcs2->crder) ||
 
1086
        !wcsutil_Eq(naxis, wcs1->csyer, wcs2->csyer) ||
 
1087
        strncmp(wcs1->dateavg, wcs2->dateavg, 72) ||
 
1088
        strncmp(wcs1->dateobs, wcs2->dateobs, 72) ||
 
1089
        wcs1->equinox != wcs2->equinox ||
 
1090
        wcs1->mjdavg != wcs2->mjdavg ||
 
1091
        wcs1->mjdobs != wcs2->mjdobs ||
 
1092
        !wcsutil_Eq(3, wcs1->obsgeo, wcs2->obsgeo) ||
 
1093
        strncmp(wcs1->radesys, wcs2->radesys, 72) ||
 
1094
        strncmp(wcs1->specsys, wcs2->specsys, 72) ||
 
1095
        strncmp(wcs1->ssysobs, wcs2->ssysobs, 72) ||
 
1096
        wcs1->velosys != wcs2->velosys ||
 
1097
        wcs1->zsource != wcs2->zsource ||
 
1098
        strncmp(wcs1->ssyssrc, wcs2->ssyssrc, 72) ||
 
1099
        wcs1->velangl != wcs2->velangl ||
 
1100
        strncmp(wcs1->wcsname, wcs2->wcsname, 72)) {
 
1101
      return 0;
 
1102
    }
 
1103
  }
 
1104
 
 
1105
  /* Compare tabular parameters */
 
1106
  if (wcs1->ntab != wcs2->ntab) {
 
1107
    return 0;
 
1108
  }
 
1109
 
 
1110
  for (i = 0; i < wcs1->ntab; ++i) {
 
1111
    if ((status = tabcmp(0, &wcs1->tab[i], &wcs2->tab[i], &tab_equal))) {
 
1112
      return status;
 
1113
    }
 
1114
    if (!tab_equal) {
 
1115
      return 0;
 
1116
    }
 
1117
  }
 
1118
 
 
1119
  *equal = 1;
 
1120
  return 0;
 
1121
}
 
1122
 
 
1123
/*--------------------------------------------------------------------------*/
 
1124
 
980
1125
int wcsfree(struct wcsprm *wcs)
981
1126
 
982
1127
{