     
     //1: Offset detector
     
-   
+    
     //S doesnt need to chagne
     
     

    float imin,imax,jmin,jmax,kmin,kmax;
    // for X
    if( source.x<pixel1D.x){
        imin=(am==axm)? 1             : ceil (source.x+am*ray.x);
        imax=(aM==axM)? geo.nVoxelX : floor(source.x+aM*ray.x);
    }else{
        imax=(am==axm)? geo.nVoxelX-1 : floor(source.x+am*ray.x);
        imin=(aM==axM)? 0             : ceil (source.x+aM*ray.x);
    }
    // for Y
    if( source.y<pixel1D.y){
        jmin=(am==aym)? 1             : ceil (source.y+am*ray.y);
        jmax=(aM==ayM)? geo.nVoxelY : floor(source.y+aM*ray.y);
    }else{
        jmax=(am==aym)? geo.nVoxelY-1 : floor(source.y+am*ray.y);
        jmin=(aM==ayM)? 0             : ceil (source.y+aM*ray.y);
    }
    // for Z
    if( source.z<pixel1D.z){
        kmin=(am==azm)? 1             : ceil (source.z+am*ray.z);
        kmax=(aM==azM)? geo.nVoxelZ : floor(source.z+aM*ray.z);
    }else{
        kmax=(am==azm)? geo.nVoxelZ-1 : floor(source.z+am*ray.z);
        kmin=(aM==azM)? 0             : ceil (source.z+aM*ray.z);
    }
    
    // get intersection point N1. eq(20-21) [(also eq 9-10)]
    float ax,ay,az;
    ax=(source.x<pixel1D.x)?  (imin-source.x)/(ray.x+0.000000000001) :  (imax-source.x)/(ray.x+0.000000000001);
    ay=(source.y<pixel1D.y)?  (jmin-source.y)/(ray.y+0.000000000001) :  (jmax-source.y)/(ray.y+0.000000000001);
    az=(source.z<pixel1D.z)?  (kmin-source.z)/(ray.z+0.000000000001) :  (kmax-source.z)/(ray.z+0.000000000001);
    
    
    
    // get index of first intersection. eq (26) and (19)
    int i,j,k;
    float aminc=min(min(ax,ay),az);
    i=(int)floor(source.x+ (aminc+am)/2*ray.x);
    j=(int)floor(source.y+ (aminc+am)/2*ray.y);
    k=(int)floor(source.z+ (aminc+am)/2*ray.z);
    // Initialize
    float ac=am;
    //eq (28), unit anlges
    float axu,ayu,azu;
    axu=1/abs(ray.x);
    ayu=1/abs(ray.y);
    azu=1/abs(ray.z);
    // eq(29), direction of update
    float iu,ju,ku;
    iu=(source.x< pixel1D.x)? 1 : -1;
    ju=(source.y< pixel1D.y)? 1 : -1;
    ku=(source.z< pixel1D.z)? 1 : -1;
    
    float maxlength=sqrt(ray.x*ray.x*geo.dVoxelX*geo.dVoxelX+ray.y*ray.y*geo.dVoxelY*geo.dVoxelY+ray.z*ray.z*geo.dVoxelZ*geo.dVoxelZ);
    float sum=0;
    unsigned int Np=(imax-imin+1)+(jmax-jmin+1)+(kmax-kmin+1); // Number of intersections
    // Go iterating over the line, intersection by intersection. If double point, no worries, 0 will be computed
    
    for (unsigned int ii=0;ii<Np;ii++){
        if (ax==aminc){
            sum+=(ax-ac)*tex3D(tex, i+0.5, j+0.5, k+0.5);
            i=i+iu;
            ac=ax;
            ax+=axu;
        }else if(ay==aminc){
            sum+=(ay-ac)*tex3D(tex, i+0.5, j+0.5, k+0.5);
            j=j+ju;
            ac=ay;
            ay+=ayu;
        }else if(az==aminc){
            sum+=(az-ac)*tex3D(tex, i+0.5, j+0.5, k+0.5);
            k=k+ku;
            ac=az;
            az+=azu;
        }
        aminc=min(min(ax,ay),az);
    }
    detector[idx]=sum*maxlength;
}


int siddon_ray_projection(float const * const img, Geometry geo, float** result,float const * const angles,int nangles){
    
    
    
    //DONE, Tesla found
    
    // copy data to CUDA memory
    cudaArray *d_imagedata = 0;
    
    const cudaExtent extent = make_cudaExtent(geo.nVoxelX, geo.nVoxelY, geo.nVoxelZ);
    cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
    cudaMalloc3DArray(&d_imagedata, &channelDesc, extent);
    if(cudaCheckErrors("cudaMalloc3D error 3D tex")){return 1;}
    
    cudaMemcpy3DParms copyParams = { 0 };
    copyParams.srcPtr = make_cudaPitchedPtr((void*)img, extent.width*sizeof(float), extent.width, extent.height);
    copyParams.dstArray = d_imagedata;
    copyParams.extent = extent;
    copyParams.kind = cudaMemcpyHostToDevice;
    cudaMemcpy3D(&copyParams);
    
    if(cudaCheckErrors("cudaMemcpy3D fail")){return 1;}
    
    // Configure texture options
    tex.normalized = false;
    tex.filterMode = cudaFilterModePoint; //we dont want interpolation
    tex.addressMode[0] = cudaAddressModeBorder;
    tex.addressMode[1] = cudaAddressModeBorder;
    tex.addressMode[2] = cudaAddressModeBorder;
    
    cudaBindTextureToArray(tex, d_imagedata, channelDesc);
    
    if(cudaCheckErrors("3D texture memory bind fail")){return 1;}
    
    
    
    
    //Done! Image put into texture memory.
    
    
    size_t num_bytes = geo.nDetecU*geo.nDetecV * sizeof(float);
    float* dProjection;
    cudaMalloc((void**)&dProjection, num_bytes);
    cudaMemset(dProjection,0,num_bytes);
    if(cudaCheckErrors("cudaMalloc fail")){return 1;}
    
    
    bool timekernel=false; // For debuggin purposes
    cudaEvent_t start, stop;
    float elapsedTime;

    Point3D source, deltaU, deltaV, uvOrigin;
    
    // 16x16 gave the best performance empirically
    // Funnily that makes it compatible with most GPUs.....
    int divU,divV;
    divU=16;
    divV=16;
    dim3 grid((geo.nDetecU+divU-1)/divU,(geo.nDetecV+divV-1)/divV,1);
    dim3 block(divU,divV,1);
    
    for (unsigned int i=0;i<nangles;i++){
        geo.alpha=angles[i*3];
        geo.theta=angles[i*3+1];
        geo.psi  =angles[i*3+2];
        //precomute distances for faster execution
        //Precompute per angle constant stuff for speed
        computeDeltas_Siddon(geo,i, &uvOrigin, &deltaU, &deltaV, &source);
        //Ray tracing!
        if (timekernel){
            cudaEventCreate(&start);
            cudaEventRecord(start,0);
        }
        kernelPixelDetector<<<grid,block>>>(geo,dProjection, source, deltaU, deltaV, uvOrigin);
        if(cudaCheckErrors("Kernel fail")){return 1;}
        if (timekernel){
            cudaEventCreate(&stop);
            cudaEventRecord(stop,0);
            cudaEventSynchronize(stop);
            cudaEventElapsedTime(&elapsedTime, start,stop);
            //mexPrintf("%f\n" ,elapsedTime);
        }
        // copy result to host
        cudaMemcpy(result[i], dProjection, num_bytes, cudaMemcpyDeviceToHost);
        if(cudaCheckErrors("cudaMemcpy fail")){return 1;}
        
        
    }
    
    
    cudaUnbindTexture(tex);
    if(cudaCheckErrors("Unbind  fail")){return 1;}
    cudaFree(dProjection);
    cudaFreeArray(d_imagedata);
    if(cudaCheckErrors("cudaFree d_imagedata fail")){return 1;}
    
    
    
    cudaDeviceReset();
    return 0;
}



/* This code precomputes The location of the source and the Delta U and delta V (in the warped space)
 * to compute the locations of the x-rays. While it seems verbose and overly-optimized,

#include <cuda_runtime_api.h>
#include <cuda.h>
#include "Siddon_projection.hpp"
#include "errors.hpp"
#include <stdio.h>
#include <math.h>

inline int cudaCheckErrors(const char * msg)
{
   cudaError_t __err = cudaGetLastError();
   if (__err != cudaSuccess)
   {
      printf("CUDA:siddo_projection:%s:%s\n",msg, cudaGetErrorString(__err));
      cudaDeviceReset();
      return 1;
   }
   return 0;
}
    
#define MAXTREADS 1024
#define PROJ_PER_BLOCK 9
#define PIXEL_SIZE_BLOCK 9
    /*GEOMETRY DEFINITION
     *
     *                Detector plane, behind
     *            |-----------------------------|
     *            |                             |
     *            |                             |
     *            |                             |
     *            |                             |
     *            |      +--------+             |
     *            |     /        /|             |
     *   A Z      |    /        / |*D           |
     *   |        |   +--------+  |             |
     *   |        |   |        |  |             |
     *   |        |   |     *O |  +             |
     *    --->y   |   |        | /              |
     *  /         |   |        |/               |
     * V X        |   +--------+                |
     *            |-----------------------------|
     *
     *           *S
     *
     *
     *
     *
     *
     **/
    
    void CreateTexture(int num_devices,const float* imagedata,Geometry geo,cudaArray** d_cuArrTex, cudaTextureObject_t *texImage,bool alloc);

__constant__ Point3D projParamsArrayDev[4*PROJ_PER_BLOCK];  // Dev means it is on device


__global__ void vecAddInPlace(float *a, float *b, unsigned long  n)
{
    int idx = blockIdx.x*blockDim.x+threadIdx.x;
    // Make sure we do not go out of bounds
    if (idx < n)
        a[idx] = a[idx] + b[idx];
}

__global__ void kernelPixelDetector( Geometry geo,
        float* detector,
        const int currProjSetNumber,
        const int totalNoOfProjections,
        cudaTextureObject_t tex){
    
    
    unsigned long y = blockIdx.y * blockDim.y + threadIdx.y;
    unsigned long x = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned long projNumber=threadIdx.z;
    
    
    if ((x>= geo.nDetecU) | (y>= geo.nDetecV)|  (projNumber>=PROJ_PER_BLOCK))
        return;
    
    size_t idx =  (size_t)(x  * geo.nDetecV + y)+ (size_t)projNumber*geo.nDetecV *geo.nDetecU ;
    int indAlpha = currProjSetNumber*PROJ_PER_BLOCK+projNumber;  // This is the ABSOLUTE projection number in the projection array

    if(indAlpha>=totalNoOfProjections)
        return;
    
    Point3D uvOrigin = projParamsArrayDev[4*projNumber];  // 6*projNumber because we have 6 Point3D values per projection
    Point3D deltaU = projParamsArrayDev[4*projNumber+1];
    Point3D deltaV = projParamsArrayDev[4*projNumber+2];
    Point3D source = projParamsArrayDev[4*projNumber+3];
    
    /////// Get coordinates XYZ of pixel UV
    int pixelV = geo.nDetecV-y-1;
