source: trunk/MESOSCALE/LMD_MM_MARS/SRC/WPS/geogrid/util/retile-cat.c @ 1243

Last change on this file since 1243 was 11, checked in by aslmd, 14 years ago

spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 16.2 KB
Line 
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
27int **** data_cache = NULL;
28char ** fname_cache = NULL;
29int * lru = NULL;
30
31int *** supertile = NULL;
32int supertile_min_x = 999999;
33int supertile_min_y = 999999;
34
35void 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
50int *** 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
96int *** 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
174void 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
462void 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
509int 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
528int 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}
Note: See TracBrowser for help on using the repository browser.