viewshed.c (e9d989e) - plDarkSky · Pawel Netzel

viewshed.c

#include <stdlib.h>
#include <math.h>
#include <omp.h>
#include <ezgdal.h>

void checkLine(EZGDAL_FRAME *frame, double **dist, int **visible, int xc, int yc, double h0, double scaleH, double r, int x1, int y1)
{
   double min = 99999999999.0;
   double a,d;
   double **buf = frame->buffer;
   EZGDAL_LAYER *l = frame->owner.stripe->layer;
   
   int dx =  abs(x1-xc), sx = xc<x1 ? 1 : -1;
   int dy = -abs(y1-yc), sy = yc<y1 ? 1 : -1;
   int err = dx+dy, e2;                                   /* error value e_xy */

if(ezgdal_is_null(l,buf[xc][yc])) return;

visible[xc][yc] = 1;
   h0 += scaleH*buf[xc][xc];

for (;;){
   
      d = dist[xc][yc];

if(d>r) break;

if(visible[xc][yc]==-1) {
        if(!ezgdal_is_null(l,buf[xc][yc])) {
          a=(h0-scaleH*buf[xc][yc])/d;
          if(a<=min) {
            visible[xc][yc] = 1;
            min=a;
          } else {
            visible[xc][yc] = 0;
          }
        }
      }
      
      e2 = 2*err;
      if (e2 >= dy) {                                         /* e_xy+e_x > 0 */
        if (xc == x1) break;
        err += dy; xc += sx;
      }
      if (e2 <= dx) {                                         /* e_xy+e_y < 0 */
         if (yc == y1) break;
         err += dx; yc += sy;
      }
   }
}

void init_mask(int **mask, int size) {
  int x, y;
#pragma omp parallel for private(x) shared(size,mask)
  for(y=0; y<size; y++) 
    for(x=0; x<size; x++)
      mask[y][x]=-1;
}

void calc_visibility_mask(EZGDAL_FRAME *height, double **distance, int **visible, double height0, double scaleH, double radius, int size) {
  int xc = (size-1)/2;
  int yc = xc;
  int x, y;

init_mask(visible,size);

#pragma omp parallel sections shared(height,distance,visible,xc,yc,height0,scaleH,radius,size)
  {
    #pragma omp section
    {
      checkLine(height,distance,visible,xc,yc,height0,scaleH,radius,0,0);
    }
    
    #pragma omp section
    {
      checkLine(height,distance,visible,xc,yc,height0,scaleH,radius,0,size-1);
    }
    
    #pragma omp section
    {
      checkLine(height,distance,visible,xc,yc,height0,scaleH,radius,size-1,size-1);
    }

#pragma omp section
    {
      checkLine(height,distance,visible,xc,yc,height0,scaleH,radius,size-1,0);
    }
  }

#pragma omp parallel sections shared(height,distance,visible,xc,yc,height0,scaleH,radius,size) private(x,y)
  {
    #pragma omp section
    {
      for(y=1; y<size-1; y++)
        checkLine(height,distance,visible,xc,yc,height0,scaleH,radius,0,y);
    }

#pragma omp section
    {
      for(x=1; x<size-1; x++)
        checkLine(height,distance,visible,xc,yc,height0,scaleH,radius,x,size-1);
    }

#pragma omp section
    {
      for(y=size-2; y>0; y--)
        checkLine(height,distance,visible,xc,yc,height0,scaleH,radius,size-1,y);
    }

#pragma omp section
    {
      for(x=size-2; x>0; x--)
        checkLine(height,distance,visible,xc,yc,height0,scaleH,radius,x,0);
    }
  }

}

void set_all_visible_mask(int **mask, int size, double **distance, double radius) {
  int x, y;

#pragma omp parallel for private(x) shared(size,mask,distance,radius)
  for(y=0; y<size; y++) 
    for(x=0; x<size; x++)
      if(distance[y][x]<=radius)
        mask[y][x]=1;
      else
        mask[y][x]=-1;
}

int **create_mask(int size, double **distance, double radius) {
  int i;
  int **d;
  
  d = malloc(size*sizeof(int *));
  for(i=0; i<size; i++)
    d[i] = malloc(size*sizeof(int));

set_all_visible_mask(d, size, distance, radius);

return d;
}

double **create_dist(int size, double scale) {
  int xc = (size-1)/2;
  int yc = xc;
  int x, y;
  double **d;
  
  d = malloc(size*sizeof(double *));
  for(y=0; y<size; y++)
    d[y] = malloc(size*sizeof(double));
    
  for(y=0; y<size; y++) 
    for(x=0; x<size; x++) 
      d[y][x]=sqrt(pow(scale*(x-xc),2)+pow(scale*(y-yc),2));
  
  return d;
}

void free_array(void **arr, int size) {
  int i;
  if(arr==NULL) return;
  
  for(i=0; i<size; i++)
    free(arr[i]);
  free(arr);
}

plDarkSky

Initial commit