geogrid/util/retile-cat.c

   1 #include <stdio.h>
   2 #include <sys/types.h>
   3 #include <sys/stat.h>
   4 #include <fcntl.h>
   5 #include <unistd.h>
   6 #include <math.h>
   7 #include <string.h>
   8
   9 #define MSG_FLAG 0x80000001
  10
  11 #define N_BYTES 1
  12 #define IN_TILE_DEGREES_LON 10.0
  13 #define IN_TILE_DEGREES_LAT 10.0
  14 #define IN_TILE_PTS_X 1200
  15 #define IN_TILE_PTS_Y 1200
  16 #define OUT_TILE_DEGREES_LON 20.0
  17 #define OUT_TILE_DEGREES_LAT 20.0
  18 #define OUT_TILE_PTS_X 600
  19 #define OUT_TILE_PTS_Y 600
  20 #define HALO_WIDTH 0
  21 #define ZDIM 1
  22 #define NCATS 24
  23 #define CATMIN 1
  24
  25 #define CACHESIZE 12
  26
  27 int **** data_cache = NULL;
  28 char ** fname_cache = NULL;
  29 int * lru = NULL;
  30
  31 int *** supertile = NULL;
  32 int supertile_min_x = 999999;
  33 int supertile_min_y = 999999;
  34
  35 void free_newtile(int *** x)
  36 {
  37   int i, j, k;
  38
  39   for(i=0; i<IN_TILE_PTS_X; i++)
  40   {
  41     for(j=0; j<IN_TILE_PTS_Y; j++)
  42     {
  43       free(x[i][j]);
  44     }
  45     free(x[i]);
  46   }
  47   free(x);
  48 }
  49
  50 int *** read_tile(char * fname)
  51 {
  52   int *** retval;
  53   int i, j, k, b;
  54   unsigned char buf[IN_TILE_PTS_X*IN_TILE_PTS_Y*ZDIM*N_BYTES];
  55   int fd;
  56
  57   retval = (int ***)malloc(sizeof(int **)*IN_TILE_PTS_X);
  58   for(i=0; i<IN_TILE_PTS_X; i++)
  59   {
  60     retval[i] = (int **)malloc(sizeof(int *)*IN_TILE_PTS_Y);
  61     for(j=0; j<IN_TILE_PTS_Y; j++)
  62     {
  63       retval[i][j] = (int *)malloc(sizeof(int)*ZDIM);
  64     }
  65   }
  66
  67   if ((fd = open(fname,O_RDONLY)) == -1)
  68   {
  69     fprintf(stderr,"Error opening source file %s\n", fname);
  70     return 0;
  71   }
  72
  73   read(fd, (void *)&buf, IN_TILE_PTS_X*IN_TILE_PTS_Y*ZDIM*N_BYTES);
  74
  75   close(fd);
  76
  77   /* Put buf into retval */
  78   for(i=0; i<IN_TILE_PTS_X; i++)
  79   {
  80     for(j=0; j<IN_TILE_PTS_Y; j++)
  81     {
  82       for(k=0; k<ZDIM; k++)
  83       {
  84         retval[i][j][k] = 0;
  85         for(b=0; b<N_BYTES; b++)
  86         {
  87           retval[i][j][k] |= buf[k*N_BYTES*IN_TILE_PTS_X*IN_TILE_PTS_Y+j*N_BYTES*IN_TILE_PTS_X+i*N_BYTES+b] << 8*(N_BYTES-b-1);
  88         }
  89       }
  90     }
  91   }
  92
  93   return retval;
  94 }
  95
  96 int *** get_tile_from_cache(int i, int j)
  97 {
  98   int ii, jj, kk, k, least, least_idx;
  99   int *** retval, *** localptr;
 100   char * fname;
 101
 102   fname = (char *)malloc(256);
 103
 104   i = (i/IN_TILE_PTS_X)*IN_TILE_PTS_X+1;
 105   j = (j/IN_TILE_PTS_Y)*IN_TILE_PTS_Y+1;
 106
 107   snprintf(fname,256,"%5.5i-%5.5i.%5.5i-%5.5i",i,i+IN_TILE_PTS_X-1,j,j+IN_TILE_PTS_Y-1);
 108
 109   /* Find out whether tile containing (i,j) is in cache */
 110   if (data_cache != NULL)
 111   {
 112     for(k=0; k<CACHESIZE; k++)
 113     {
 114       if (fname_cache[k] != NULL)
 115       {
 116         if (strncmp(fname_cache[k],fname,256) == 0)
 117         {
 118           free(fname);
 119           retval = data_cache[k];
 120           return retval;
 121         }
 122       }
 123     }
 124   }
 125
 126   /* If not, read from file */
 127   localptr = read_tile(fname);
 128
 129   /* Also store tile in the cache */
 130   if (data_cache == NULL)
 131   {
 132     data_cache = (int ****)malloc(sizeof(int ***)*CACHESIZE);
 133     fname_cache = (char **)malloc(sizeof(char *)*CACHESIZE);
 134     lru = (int *)malloc(sizeof(int)*CACHESIZE);
 135     for(k=0; k<CACHESIZE; k++)
 136     {
 137       data_cache[k] = NULL;
 138       fname_cache[k] = NULL;
 139       lru[k] = 0;
 140     }
 141   }
 142
 143   least = 0;
 144   least_idx = 0;
 145   for(k=0; k<CACHESIZE; k++)
 146   {
 147     lru[k]++;
 148     if (lru[k] > least)
 149     {
 150       least = lru[k];
 151       least_idx = k;
 152     }
 153   }
 154
 155   if (data_cache[least_idx] == NULL)
 156   {
 157     data_cache[least_idx] = localptr;
 158     fname_cache[least_idx] = fname;
 159     lru[least_idx] = 0;
 160   }
 161   else
 162   {
 163     free_newtile(data_cache[least_idx]);
 164     data_cache[least_idx] = localptr;
 165     free(fname_cache[least_idx]);
 166     fname_cache[least_idx] = fname;
 167     lru[least_idx] = 0;
 168   }
 169
 170   retval = localptr;
 171   return retval;
 172 }
 173
 174 void build_supertile(int i, int j)
 175 {
 176   int ii, jj, kk;
 177   int doflip;
 178   int *** newtile;
 179
 180   if (i < 0)
 181     i = i + IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
 182   if (i >= IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON))
 183     i = i - IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
 184   if (j < 0)
 185     j = j + IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
 186   if (j >= IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT))
 187     j = j - IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
 188
 189   if (supertile == NULL)
 190   {
 191     supertile = (int ***)malloc(sizeof(int **)*3*IN_TILE_PTS_X);
 192     for(ii=0; ii<3*IN_TILE_PTS_X; ii++)
 193     {
 194       supertile[ii] = (int **)malloc(sizeof(int *)*3*IN_TILE_PTS_Y);
 195       for(jj=0; jj<3*IN_TILE_PTS_Y; jj++)
 196       {
 197         supertile[ii][jj] = (int *)malloc(sizeof(int)*ZDIM);
 198       }
 199     }
 200   }
 201
 202   supertile_min_x = (i / IN_TILE_PTS_X)*IN_TILE_PTS_X;
 203   supertile_min_y = (j / IN_TILE_PTS_Y)*IN_TILE_PTS_Y;
 204
 205   /* Get tile containing (i,j) from cache*/
 206   /* Get surrounding tiles from cache */
 207
 208   /* Lower-left */
 209   ii = i - IN_TILE_PTS_X;
 210   if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 211   doflip = 0;
 212   jj = j - IN_TILE_PTS_Y;
 213   if (jj < 0)
 214   {
 215     doflip = 1;
 216     jj = j;
 217     ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 218     if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 219   }
 220   newtile = get_tile_from_cache(ii, jj);
 221   if (doflip)
 222   {
 223     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 224     {
 225       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 226       {
 227         for(kk=0; kk<ZDIM; kk++)
 228           supertile[ii][jj][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
 229       }
 230     }
 231   }
 232   else
 233   {
 234     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 235     {
 236       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 237       {
 238         for(kk=0; kk<ZDIM; kk++)
 239           supertile[ii][jj][kk] = newtile[ii][jj][kk];
 240       }
 241     }
 242   }
 243
 244   /* Left */
 245   ii = i - IN_TILE_PTS_X;
 246   if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 247   jj = j;
 248   newtile = get_tile_from_cache(ii, jj);
 249   for(ii=0; ii<IN_TILE_PTS_X; ii++)
 250   {
 251     for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 252     {
 253       for(kk=0; kk<ZDIM; kk++)
 254         supertile[ii][jj+IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
 255     }
 256   }
 257
 258   /* Upper-left */
 259   ii = i - IN_TILE_PTS_X;
 260   if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 261   doflip = 0;
 262   jj = j + IN_TILE_PTS_Y;
 263   if (jj >= (int)(180./IN_TILE_DEGREES_LAT)*IN_TILE_PTS_Y)
 264   {
 265     doflip = 1;
 266     jj = j;
 267     ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 268     if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 269   }
 270   newtile = get_tile_from_cache(ii, jj);
 271   if (doflip)
 272   {
 273     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 274     {
 275       for(jj=0; jj<IN_TILE_PTS_X; jj++)
 276       {
 277         for(kk=0; kk<ZDIM; kk++)
 278           supertile[ii][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
 279       }
 280     }
 281   }
 282   else
 283   {
 284     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 285     {
 286       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 287       {
 288         for(kk=0; kk<ZDIM; kk++)
 289           supertile[ii][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
 290       }
 291     }
 292   }
 293
 294   /* Below */
 295   ii = i;
 296   doflip = 0;
 297   jj = j - IN_TILE_PTS_Y;
 298   if (jj < 0)
 299   {
 300     doflip = 1;
 301     jj = j;
 302     ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 303     if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 304   }
 305   newtile = get_tile_from_cache(ii, jj);
 306   if (doflip)
 307   {
 308     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 309     {
 310       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 311       {
 312         for(kk=0; kk<ZDIM; kk++)
 313           supertile[ii+IN_TILE_PTS_X][jj][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
 314       }
 315     }
 316   }
 317   else
 318   {
 319     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 320     {
 321       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 322       {
 323         for(kk=0; kk<ZDIM; kk++)
 324           supertile[ii+IN_TILE_PTS_X][jj][kk] = newtile[ii][jj][kk];
 325       }
 326     }
 327   }
 328
 329   /* Center */
 330   newtile = get_tile_from_cache(i, j);
 331   for(ii=0; ii<IN_TILE_PTS_X; ii++)
 332   {
 333     for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 334     {
 335       for(kk=0; kk<ZDIM; kk++)
 336         supertile[ii+IN_TILE_PTS_X][jj+IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
 337     }
 338   }
 339
 340   /* Above */
 341   ii = i;
 342   doflip = 0;
 343   jj = j + IN_TILE_PTS_Y;
 344   if (jj >= (int)(180./IN_TILE_DEGREES_LAT)*IN_TILE_PTS_Y)
 345   {
 346     doflip = 1;
 347     jj = j;
 348     ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 349     if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 350   }
 351   newtile = get_tile_from_cache(ii, jj);
 352   if (doflip)
 353   {
 354     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 355     {
 356       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 357       {
 358         for(kk=0; kk<ZDIM; kk++)
 359           supertile[ii+IN_TILE_PTS_X][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
 360       }
 361     }
 362   }
 363   else
 364   {
 365     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 366     {
 367       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 368       {
 369         for(kk=0; kk<ZDIM; kk++)
 370           supertile[ii+IN_TILE_PTS_X][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
 371       }
 372     }
 373   }
 374
 375   /* Lower-right */
 376   ii = i + IN_TILE_PTS_X;
 377   if (ii >= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X) ii -= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 378   doflip = 0;
 379   jj = j - IN_TILE_PTS_Y;
 380   if (jj < 0)
 381   {
 382     doflip = 1;
 383     jj = j;
 384     ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 385     if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 386   }
 387   newtile = get_tile_from_cache(ii, jj);
 388   if (doflip)
 389   {
 390     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 391     {
 392       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 393       {
 394         for(kk=0; kk<ZDIM; kk++)
 395           supertile[ii+2*IN_TILE_PTS_X][jj][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
 396       }
 397     }
 398   }
 399   else
 400   {
 401     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 402     {
 403       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 404       {
 405         for(kk=0; kk<ZDIM; kk++)
 406           supertile[ii+2*IN_TILE_PTS_X][jj][kk] = newtile[ii][jj][kk];
 407       }
 408     }
 409   }
 410
 411   /* Right */
 412   ii = i + IN_TILE_PTS_X;
 413   if (ii >= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X) ii -= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 414   jj = j;
 415   newtile = get_tile_from_cache(ii, jj);
 416   for(ii=0; ii<IN_TILE_PTS_X; ii++)
 417   {
 418     for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 419     {
 420       for(kk=0; kk<ZDIM; kk++)
 421         supertile[ii+2*IN_TILE_PTS_X][jj+IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
 422     }
 423   }
 424
 425   /* Upper-right */
 426   ii = i + IN_TILE_PTS_X;
 427   if (ii >= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X) ii -= (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 428   doflip = 0;
 429   jj = j + IN_TILE_PTS_Y;
 430   if (jj >= (int)(180./IN_TILE_DEGREES_LAT)*IN_TILE_PTS_Y)
 431   {
 432     doflip = 1;
 433     jj = j;
 434     ii -= (int)(180./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 435     if (ii < 0) ii += (int)(360./IN_TILE_DEGREES_LON)*IN_TILE_PTS_X;
 436   }
 437   newtile = get_tile_from_cache(ii, jj);
 438   if (doflip)
 439   {
 440     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 441     {
 442       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 443       {
 444         for(kk=0; kk<ZDIM; kk++)
 445           supertile[ii+2*IN_TILE_PTS_X][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][IN_TILE_PTS_Y-1-jj][kk];
 446       }
 447     }
 448   }
 449   else
 450   {
 451     for(ii=0; ii<IN_TILE_PTS_X; ii++)
 452     {
 453       for(jj=0; jj<IN_TILE_PTS_Y; jj++)
 454       {
 455         for(kk=0; kk<ZDIM; kk++)
 456           supertile[ii+2*IN_TILE_PTS_X][jj+2*IN_TILE_PTS_Y][kk] = newtile[ii][jj][kk];
 457       }
 458     }
 459   }
 460 }
 461
 462 void get_value(int * num_cats, int i, int j, int k, int irad)
 463 {
 464   int i_src, j_src;
 465   int ii, jj;
 466   int n, sum;
 467   float r_rel;
 468
 469   if (i < 0)
 470     i = i + IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
 471   if (i >= IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON))
 472     i = i - IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
 473   if (j < 0)
 474     j = j + IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
 475   if (j >= IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT))
 476     j = j - IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
 477
 478   i_src = i % IN_TILE_PTS_X + IN_TILE_PTS_X;
 479   j_src = j % IN_TILE_PTS_Y + IN_TILE_PTS_Y;
 480
 481   /* Interpolate values from supertile */
 482   for(ii=0; ii<NCATS; ii++)
 483     num_cats[ii] = 0;
 484
 485   for(ii=i_src-irad; ii<=i_src+irad; ii++)
 486   {
 487     for(jj=j_src-irad; jj<=j_src+irad; jj++)
 488     {
 489       num_cats[supertile[ii][jj][k]-CATMIN]++;
 490     }
 491   }
 492 /*
 493   sum = 0;
 494   n = 0;
 495   for(ii=i_src-irad; ii<=i_src+irad; ii++)
 496   {
 497     for(jj=j_src-irad; jj<=j_src+irad; jj++)
 498     {
 499       sum += supertile[i_src][j_src][k];
 500       n++;
 501     }
 502   }
 503
 504   if (n > 0) return sum/n;
 505   else return 0;
 506 */
 507 }
 508
 509 int is_in_supertile(int i, int j)
 510 {
 511   if (i < 0)
 512     i = i + IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
 513   if (i >= IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON))
 514     i = i - IN_TILE_PTS_X*(int)(360./IN_TILE_DEGREES_LON);
 515   if (j < 0)
 516     j = j + IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
 517   if (j >= IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT))
 518     j = j - IN_TILE_PTS_Y*(int)(180./IN_TILE_DEGREES_LAT);
 519
 520   /* Check whether (i,j) is in the interior of supertile */
 521   if ((i >= supertile_min_x) && (i < supertile_min_x+IN_TILE_PTS_X) &&
 522       (j >= supertile_min_y) && (j < supertile_min_y+IN_TILE_PTS_Y))
 523     return 1;
 524   else
 525     return 0;
 526 }
 527
 528 int main(int argc, char ** argv)
 529 {
 530   int tile_x, tile_y, input_x, input_y, temp_x, temp_y, temp;
 531   int i, j, k, z, ii, jj;
 532   int i_src, j_src;
 533   int *** intdata;
 534   int * num_cats;
 535   int out_fd;
 536   int ir_rel;
 537   float r_rel;
 538   unsigned char * outdata;
 539   char out_filename[256];
 540
 541   r_rel = (float)IN_TILE_PTS_X/(float)OUT_TILE_PTS_X*OUT_TILE_DEGREES_LON/IN_TILE_DEGREES_LON;
 542   ir_rel = (int)rint(r_rel);
 543
 544   /* Allocate memory to hold a single output tile */
 545   intdata = (int ***)malloc(sizeof(int **)*(OUT_TILE_PTS_X+2*HALO_WIDTH));
 546   for(i=0; i<(OUT_TILE_PTS_X+2*HALO_WIDTH); i++)
 547   {
 548     intdata[i] = (int **)malloc(sizeof(int *)*(OUT_TILE_PTS_Y+2*HALO_WIDTH));
 549     for(j=0; j<(OUT_TILE_PTS_Y+2*HALO_WIDTH); j++)
 550     {
 551       intdata[i][j] = (int *)malloc(sizeof(int)*ZDIM);
 552     }
 553   }
 554
 555   num_cats = (int *)malloc(sizeof(int)*NCATS);
 556
 557   /* Allocate output buffer */
 558   outdata = (unsigned char *)malloc((OUT_TILE_PTS_X+2*HALO_WIDTH)*(OUT_TILE_PTS_Y+2*HALO_WIDTH)*ZDIM*N_BYTES);
 559
 560   for(tile_x=0; tile_x<OUT_TILE_PTS_X*(int)(360.0/OUT_TILE_DEGREES_LON); tile_x+=OUT_TILE_PTS_X)
 561   {
 562     for(tile_y=0; tile_y<OUT_TILE_PTS_Y*(int)(180.0/OUT_TILE_DEGREES_LAT); tile_y+=OUT_TILE_PTS_Y)
 563     {
 564       /* Build name of output file for current tile_x and tile_y */
 565       sprintf(out_filename,"retiled/%5.5i-%5.5i.%5.5i-%5.5i",tile_x+1,tile_x+OUT_TILE_PTS_X,tile_y+1,tile_y+OUT_TILE_PTS_Y);
 566
 567       /* Initialize the output data for current tile */
 568       for(i=0; i<(OUT_TILE_PTS_X+2*HALO_WIDTH); i++)
 569       {
 570         for(j=0; j<(OUT_TILE_PTS_Y+2*HALO_WIDTH); j++)
 571         {
 572             intdata[i][j][0] = MSG_FLAG;
 573         }
 574       }
 575
 576       /* Attempt to open output file */
 577       if ((out_fd = open(out_filename, O_WRONLY|O_CREAT|O_TRUNC, S_IRUSR|S_IWUSR|S_IRGRP|S_IROTH)) == -1)
 578       {
 579         fprintf(stderr, "Error: Could not create or open file %s\n", out_filename);
 580         return 1;
 581       }
 582
 583       /* Fill tile with data */
 584       for(i=-1*HALO_WIDTH; i<OUT_TILE_PTS_X+HALO_WIDTH; i++)
 585       {
 586         for(j=-1*HALO_WIDTH; j<OUT_TILE_PTS_Y+HALO_WIDTH; j++)
 587         {
 588           if (intdata[i+HALO_WIDTH][j+HALO_WIDTH][0] == MSG_FLAG)
 589           {
 590             i_src = ir_rel*(tile_x+i);
 591             j_src = ir_rel*(tile_y+j);
 592
 593             build_supertile(i_src,j_src);
 594
 595             for(ii=-1*HALO_WIDTH; ii<OUT_TILE_PTS_X+HALO_WIDTH; ii++)
 596             {
 597               for(jj=-1*HALO_WIDTH; jj<OUT_TILE_PTS_Y+HALO_WIDTH; jj++)
 598               {
 599                 i_src = ir_rel*(tile_x+ii);
 600                 j_src = ir_rel*(tile_y+jj);
 601                 if (is_in_supertile(i_src,j_src))
 602                 {
 603                   get_value(num_cats, i_src, j_src, 0, ir_rel-1);
 604                   intdata[ii+HALO_WIDTH][jj+HALO_WIDTH][0] = CATMIN-1;
 605                   temp = 0;
 606                   for(k=0; k<NCATS; k++)
 607                   {
 608                     if (num_cats[k] >= temp)
 609                     {
 610                       temp = num_cats[k];
 611                       intdata[ii+HALO_WIDTH][jj+HALO_WIDTH][0] = k+CATMIN;
 612                     }
 613                   }
 614                 }
 615               }
 616             }
 617
 618           }
 619         }
 620       }
 621
 622       /* Write out the data */
 623       for(i=0; i<(OUT_TILE_PTS_X+2*HALO_WIDTH); i++)
 624       {
 625         for(j=0; j<(OUT_TILE_PTS_Y+2*HALO_WIDTH); j++)
 626         {
 627           for(z=0; z<ZDIM; z++)
 628           {
 629             for(k=0; k<N_BYTES; k++)
 630             {
 631               outdata[z*N_BYTES*(OUT_TILE_PTS_Y+2*HALO_WIDTH)*(OUT_TILE_PTS_X+2*HALO_WIDTH)+j*N_BYTES*(OUT_TILE_PTS_X+2*HALO_WIDTH)+i*N_BYTES+k] =
 632                      (intdata[i][j][z] & (0xff << 8*(N_BYTES-k-1))) >> (8*(N_BYTES-k-1));
 633             }
 634           }
 635         }
 636       }
 637       write(out_fd,(void *)outdata,(OUT_TILE_PTS_X+2*HALO_WIDTH)*(OUT_TILE_PTS_Y+2*HALO_WIDTH)*ZDIM*N_BYTES);
 638       close(out_fd);
 639       printf("Wrote file %s\n",out_filename);
 640     }
 641   }
 642
 643   free(num_cats);
 644
 645   /* Deallocate memory */
 646   for(i=0; i<(OUT_TILE_PTS_X+2*HALO_WIDTH); i++)
 647   {
 648     for(j=0; j<(OUT_TILE_PTS_Y+2*HALO_WIDTH); j++)
 649     {
 650       free(intdata[i][j]);
 651     }
 652     free(intdata[i]);
 653   }
 654   free(intdata);
 655
 656
 657   return 0;
 658 }