source: trunk/MESOSCALE/LMD_MM_MARS/SRC/WPS/geogrid/util/retile-cont.c @ 3610

Last change on this file since 3610 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 OUT_ZDIM 24
23#define NCATS 24
24#define CATMIN 1
25
26#define CACHESIZE 12
27
28int **** data_cache = NULL;
29char ** fname_cache = NULL;
30int * lru = NULL;
31
32int *** supertile = NULL;
33int supertile_min_x = 999999;
34int supertile_min_y = 999999;
35
36void 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
51int *** 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
97int *** 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
175void 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
463void 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
510int 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
529int 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}
Note: See TracBrowser for help on using the repository browser.