[postgis-users] Adding a arbitrary ray function to Postgis raster

[Nathaniel Hunter Clay]

I am writing a function to add the ability get a ray of pixels from a 
raster. To facilitate functions such as line of sight, viewshed and 
other radio wave propagation analysis. The main reason for writing the 
function in c is to ensure that, every pixel along the ray is returned 
and for the use in a future addition of a call back function that 
returns every ray from a centre point to the edge of the raster. However 
my c is a little rusty and I am getting a segmentation fault.

I have a function that I am passing a pointer to  then allocating memory 
inside the function.

     count = rt_band_get_ray(band,x1,y1,x2,y2,&npixels);

Here is the function body:

     /**
     * Processes an array of pixels that fall on the ray from x1,y1 to 
x2,y2.
      * @param band : band to be processed
      * @param x1 : starting x-cordinate ( 0 based )
      * @param y1 : starting y-cordinate ( 0 based )
      * @param x2 : ending x-crodinate ( 0 based )
      * @param y2 : ending y-cordinate ( 0 based )
      * @param npixels : pointer to array of pixels
      *
      * @return -1 on error 0 for success.
      * http://en.wikipedia.org/wiki/Bresenham%27s_line_algorithm
      */
     int
     rt_band_get_ray(rt_band band, int x1, int y1, int x2, int y2, 
rt_pixel *npixels )
     {
       int x;
       int y;
       int xinc1;
       int yinc1;
       int xinc2;
       int yinc2;
       int dx;
       int dy;
       int den;
       int num;
       int numadd;
       int curpixel;
       double pixval;
       int isnodata = 0;
       rt_pixel pixel = NULL;
       int count;

       dx = abs( x2 - x1 );
       dy = abs( y2 - y1 );

       x = x1;
       y = y1;

       // considerations for all eight quads

       if ( x2 >= x1 )  // the x-values are increasing
         {
           xinc1 = 1;
           xinc2 = 1;
         }
       else   // the x-values are decreasing
         {
         xinc1 = -1;
         xinc2 = -1;
         }

       if( y2 >= y1 ) // the y-values are increasing
         {
           yinc1 = 1;
           yinc2 = 1;
         }
       else  // the y-values are decreasing
         {
           yinc1 = -1;
           yinc2 = -1;
         }

       if ( dx >= dy )  // there is at least one x-value for every y-value
         {
           xinc1 = 0;   //don't change x when numerator >= denominator
           yinc2 = 0;   // don't change y for every iteration
           den = dx;
           num = dx / 2;
           numadd = dy;
           count = dx;  // there are more x-values than y-values
         }
       else
         {
           xinc2 = 0;
           yinc1 = 0;
           den = dy;
           num = dy / 2;
           numadd = dx;
           count = dy;
         }

       // allocate array
         if (npixels == NULL)
          *npixels = (rt_pixel) rtalloc(sizeof(struct rt_pixel_t) * count);
         else
          *npixels = (rt_pixel) rtrealloc(*npixels, sizeof(struct 
rt_pixel_t) * count);
         if (npixels == NULL) {
           rterror("rt_band_get_ray: Unable to allocate memory for ray 
pixel(s)");
           return -1;
         }

       for ( curpixel = 0; curpixel < count; curpixel++) //itterate 
through all pixels on the line
         {
           // set current pixel to array
           rt_band_get_pixel(band,x,y,&pixval,&isnodata);
           pixel = &((*npixels)[curpixel]);
           pixel->x = x;
           pixel->y = y;
           pixel->nodata = isnodata;
           pixel->value = pixval;

           num += numadd;  // increase the numerator by the top of the 
fraction
           if ( num >= den )  // check if numerator >= denominator
         {
           num -= den;  // calculate the new numerator value
           x += xinc1;  // change x as appropriate
           y += yinc1;  // change y as appropriate
         }
           x += xinc2;   // change x as appropriate
           y += yinc2;   // change y as appropriate
         }
       return count; // return number of npixel elements.

     }
I think I am loosing scope on the memory allocated inside the function. 
And experiencing a segmentation fault. On this line:


     for ( i = 0; i < count; i++ ){
-->          pixel = &((*npixels)[i]);
           results[i] =  pixel->value;
         }

Am I loosing scope, on the data allocated inside the function? If so how 
do I accomplish, allocating an array inside this function with out 
loosing scope?

I humbly submit to any questions, comments or suggestions to fixing my 
code or to the entire idea its self.

Humbly,

Nathaniel Hunter Clay