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