     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);
+        imin=(am==axm)? 1.0f             : ceilf (source.x+am*ray.x);
+        imax=(aM==axM)? geo.nVoxelX      : floorf(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);
+        imax=(am==axm)? geo.nVoxelX-1.0f : floorf(source.x+am*ray.x);
+        imin=(aM==axM)? 0.0f             : ceilf (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);
+        jmin=(am==aym)? 1.0f             : ceilf (source.y+am*ray.y);
+        jmax=(aM==ayM)? geo.nVoxelY      : floorf(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);
+        jmax=(am==aym)? geo.nVoxelY-1.0f : floorf(source.y+am*ray.y);
+        jmin=(aM==ayM)? 0.0f             : ceilf (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);
+        kmin=(am==azm)? 1.0f             : ceilf (source.z+am*ray.z);
+        kmax=(aM==azM)? geo.nVoxelZ      : floorf(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);
+        kmax=(am==azm)? geo.nVoxelZ-1.0f : floorf(source.z+am*ray.z);
+        kmin=(aM==azM)? 0.0f             : ceilf (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);
+    ax=(source.x<pixel1D.x)?  __fdividef(imin-source.x,ray.x+0.000000000001f) :  __fdividef(imax-source.x,ray.x+0.000000000001f);
+    ay=(source.y<pixel1D.y)?  __fdividef(jmin-source.y,ray.y+0.000000000001f) :  __fdividef(jmax-source.y,ray.y+0.000000000001f);
+    az=(source.z<pixel1D.z)?  __fdividef(kmin-source.z,ray.z+0.000000000001f) :  __fdividef(kmax-source.z,ray.z+0.000000000001f);
     
     
     
     // 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);
+    float aminc=fminf(fminf(ax,ay),az);
+    i=(int)floorf(source.x+ (aminc+am)*0.5f*ray.x);
+    j=(int)floorf(source.y+ (aminc+am)*0.5f*ray.y);
+    k=(int)floorf(source.z+ (aminc+am)*0.5f*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);
+    axu=__frcp_rd(fabsf(ray.x));
+    ayu=__frcp_rd(fabsf(ray.y));
+    azu=__frcp_rd(fabsf(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;
+    iu=(source.x< pixel1D.x)? 1.0f : -1.0f;
+    ju=(source.y< pixel1D.y)? 1.0f : -1.0f;
+    ku=(source.z< pixel1D.z)? 1.0f : -1.0f;
+    
+    float maxlength=__fsqrt_rd(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.0f;
     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
-    
+    i+=0.5f;
+    j+=0.5f;
+    k+=0.5f;
     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);
+            sum+=(ax-ac)*tex3D<float>(tex, i, j, k);
             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);
+            sum+=(ay-ac)*tex3D<float>(tex, i, j, k);
             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);
+            sum+=(az-ac)*tex3D<float>(tex, i, j, k);
             k=k+ku;
             ac=az;
             az+=azu;
         }
-        aminc=min(min(ax,ay),az);
+        aminc=fminf(fminf(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;
-    
+int siddon_ray_projection(float  *  img, Geometry geo, float** result,float const * const angles,int nangles){
+
+    
+    // Prepare for MultiGPU
+    int deviceCount = 0;
+    cudaGetDeviceCount(&deviceCount);
+    if(cudaCheckErrors("Device query fail")){return 1;}
+    if (deviceCount == 0) {
+        return ERR_NO_CAPABLE_DEVICES;
+        }
+    //
+    // CODE assumes
+    // 1.-All available devices are usable by this code
+    // 2.-All available devices are equal, they are the same machine (warning thrown)
+    int dev;
+    char * devicenames;
+    cudaDeviceProp deviceProp;
+    
+    for (dev = 0; dev < deviceCount; dev++) {
+        cudaSetDevice(dev);
+        cudaGetDeviceProperties(&deviceProp, dev);
+        if (dev>0){
+            if (strcmp(devicenames,deviceProp.name)!=0){
+                printf("Ax:GPUselect","Detected one (or more) different GPUs.\n This code is not smart enough to separate the memory GPU wise if they have different computational times or memory limits.\n First GPU parameters used. If the code errors you might need to change the way GPU selection is performed. \n Siddon_projection.cu line 275.");
+                break;
+            }
+        }
+        devicenames=deviceProp.name;
+    }
+    
+    
+    
+    
+    // Check free memory
+    size_t mem_GPU_global;
+    checkFreeMemory(deviceCount,&mem_GPU_global);
+    if(cudaCheckErrors("Check mem error")){return ERR_NO_FREE_DEVICES;}
+
+    size_t mem_image=                 (unsigned long long)geo.nVoxelX*(unsigned long long)geo.nVoxelY*(unsigned long long)geo.nVoxelZ*sizeof(float);
+    size_t mem_proj=                  (unsigned long long)geo.nDetecU*(unsigned long long)geo.nDetecV*sizeof(float);
+    
+    // Does everything fit in the GPUs?
+    bool fits_in_memory=false;
+    unsigned int splits=1;
+    Geometry * geoArray;
+    
+    
+    if (mem_image+2*PROJ_PER_BLOCK*mem_proj<mem_GPU_global){// yes it does
+        fits_in_memory=true;
+        geoArray=(Geometry*)malloc(sizeof(Geometry));
+        geoArray[0]=geo;
+    }
+    else{// Nope nope.
+        fits_in_memory=false; // Oh dear.
+        // approx free memory we have. We already have left some extra 5% free for internal stuff
+        // we need a second projection memory to combine multi-GPU stuff.
+        size_t mem_free=mem_GPU_global-4*PROJ_PER_BLOCK*mem_proj;
+        
+        
+        splits=mem_image/mem_free+1;// Ceil of the truncation
+        geoArray=(Geometry*)malloc(splits*sizeof(Geometry));
+        splitImage(splits,geo,geoArray,nangles);
+    }
+    
+    // Allocate axuiliary memory for projections on the GPU to accumulate partial results
+    float ** dProjection_accum;
+    size_t num_bytes_proj = PROJ_PER_BLOCK*geo.nDetecU*geo.nDetecV * sizeof(float);
+    if (!fits_in_memory){
+        dProjection_accum=(float**)malloc(2*deviceCount*sizeof(float*));
+        for (dev = 0; dev < deviceCount; dev++) {
+            cudaSetDevice(dev);
+            for (int i = 0; i < 2; ++i){
+                cudaMalloc((void**)&dProjection_accum[dev*deviceCount+i], num_bytes_proj);
+                cudaMemset(dProjection_accum[dev*deviceCount+i],0,num_bytes_proj);
+                if(cudaCheckErrors("cudaMallocauxiliarty projections fail")){return 1;}
+            }
+        }
+    }
+    
+    // This is happening regarthless if the image fits on memory
+    float** dProjection=(float**)malloc(2*deviceCount*sizeof(float*));
+    for (dev = 0; dev < deviceCount; dev++){
+        cudaSetDevice(dev);
+        
+        for (int i = 0; i < 2; ++i){
+            cudaMalloc((void**)&dProjection[dev*deviceCount+i],   num_bytes_proj);
+            cudaMemset(dProjection[dev*deviceCount+i]  ,0,num_bytes_proj);
+            if(cudaCheckErrors("cudaMalloc projections fail")){return 1;}
+        }
+    }
+    
+    
+    //Pagelock memory for syncronous copy.
+    // Lets try to make the host memory pinned:
+    // We laredy queried the GPU and assuemd they are the same, thus shoudl have the same attributes.
+    int isHostRegisterSupported;
+    cudaDeviceGetAttribute(&isHostRegisterSupported,cudaDevAttrHostRegisterSupported,0);
+    // empirical testing shows that when the image split is smaller than 1 (also implies the image is not very big), the time to
+    // pin the memory is greater than the lost time in Syncronously launching the memcpys. This is only worth it when the image is too big.
+    if (isHostRegisterSupported & splits>1){
+        cudaHostRegister(img, (size_t)geo.nVoxelX*(size_t)geo.nVoxelY*(size_t)geo.nVoxelZ*(size_t)sizeof(float),cudaHostRegisterPortable);
+    }
+
+    
+    
+    // auxiliary variables
+    Point3D source, deltaU, deltaV, uvOrigin;
+    Point3D* projParamsArrayHost;
+    cudaMallocHost((void**)&projParamsArrayHost,4*PROJ_PER_BLOCK*sizeof(Point3D));
+    
+    
+    // Create Streams for overlapping memcopy and compute
+    int nStreams=deviceCount*2;
+    cudaStream_t* stream=(cudaStream_t*)malloc(nStreams*sizeof(cudaStream_t));;
+    
+    
+    for (dev = 0; dev < deviceCount; dev++){
+        cudaSetDevice(dev);
+        for (int i = 0; i < 2; ++i){
+            cudaStreamCreate(&stream[i+dev*2]);
+            
+        }
+    }
+    int nangles_device=(nangles+deviceCount-1)/deviceCount;
+    if(cudaCheckErrors("Stream creation fail")){return 1;}
+    
+    cudaTextureObject_t *texImg = new cudaTextureObject_t[deviceCount];
+    cudaArray **d_cuArrTex = new cudaArray*[deviceCount];
+    
+    for (unsigned int sp=0;sp<splits;sp++){
+        
+        // Create texture objects for all GPUs
+        
+        
+        size_t linear_idx_start;
+        //First one shoudl always be  the same size as all the rest but the last
+        linear_idx_start= (size_t)sp*(size_t)geoArray[0].nVoxelX*(size_t)geoArray[0].nVoxelY*(size_t)geoArray[0].nVoxelZ;
+        
+        
+        CreateTexture(deviceCount,&img[linear_idx_start],geoArray[sp],d_cuArrTex,texImg,!sp);
+        if(cudaCheckErrors("Texture object creation fail")){return 1;}
+        
+        
+        // Prepare kernel lauch variables
+        
+        int divU,divV;
+        divU=PIXEL_SIZE_BLOCK;
+        divV=PIXEL_SIZE_BLOCK;
+        dim3 grid((geoArray[sp].nDetecU+divU-1)/divU,(geoArray[0].nDetecV+divV-1)/divV,1);
+        dim3 block(divU,divV,PROJ_PER_BLOCK);
+        
+        unsigned int proj_global;
+        unsigned int noOfKernelCalls = (nangles_device+PROJ_PER_BLOCK-1)/PROJ_PER_BLOCK;  // We'll take care of bounds checking inside the loop if nalpha is not divisible by PROJ_PER_BLOCK
+        unsigned int i;
+        // Now that we have prepared the image (piece of image) and parameters for kernels
+        // we project for all angles.
+        for ( i=0; i<noOfKernelCalls; i++){
+            for (dev=0;dev<deviceCount;dev++){
+                cudaSetDevice(dev);
+                
+                for(unsigned int j=0; j<PROJ_PER_BLOCK; j++){
+                    proj_global=(i*PROJ_PER_BLOCK+j)+dev*nangles_device;
+                    if (proj_global>=nangles)
+                        break;
+                    geoArray[sp].alpha=angles[proj_global*3];
+                    geoArray[sp].theta=angles[proj_global*3+1];
+                    geoArray[sp].psi  =angles[proj_global*3+2];
+                    
+                    
+                    //precomute distances for faster execution
+                    //Precompute per angle constant stuff for speed
+                    computeDeltas_Siddon(geoArray[sp],proj_global, &uvOrigin, &deltaU, &deltaV, &source);
+                    //Ray tracing!
+                    projParamsArrayHost[4*j]=uvOrigin;		// 6*j because we have 6 Point3D values per projection
+                    projParamsArrayHost[4*j+1]=deltaU;
+                    projParamsArrayHost[4*j+2]=deltaV;
+                    projParamsArrayHost[4*j+3]=source;
+                    
+                }
+                cudaMemcpyToSymbolAsync(projParamsArrayDev, projParamsArrayHost, sizeof(Point3D)*4*PROJ_PER_BLOCK,0,cudaMemcpyHostToDevice,stream[dev*2]);
+                cudaStreamSynchronize(stream[dev*2]);
+                if(cudaCheckErrors("kernel fail")){return 1;}
+                kernelPixelDetector<<<grid,block,0,stream[dev*2]>>>(geoArray[sp],dProjection[(i%2)+dev*deviceCount],i,nangles,texImg[dev]);
+            }
+            // Now that the computation is happening, we need to either prepare the memory for
+            // combining of the projections (splits>1) or start removing previous results.
+
+            // If our image does not fit in memory then we need to make sure we accumulate previous results too.
+            if( !fits_in_memory){
+                if(sp>0){
+                    // First, grab previous results and put them in the auxiliary variable
+                    for (dev = 0; dev < deviceCount; dev++){
+                        cudaSetDevice(dev);
+                        cudaMemcpyAsync(dProjection_accum[(i%2)+dev*deviceCount], result[i*PROJ_PER_BLOCK+dev*nangles_device], num_bytes_proj, cudaMemcpyHostToDevice,stream[dev*2+1]);
+                    }            
+                    // Second, take the results from current compute call and add it to the code in execution.
+                    for (dev = 0; dev < deviceCount; dev++){
+                        cudaSetDevice(dev);
+                        cudaStreamSynchronize(stream[dev*2+1]); // wait until copy is finished
+                        vecAddInPlace<<<(geo.nDetecU*geo.nDetecV*PROJ_PER_BLOCK+MAXTREADS-1)/MAXTREADS,MAXTREADS,0,stream[dev*2]>>>(dProjection[(i%2)+dev*deviceCount],dProjection_accum[(i%2)+dev*deviceCount],(unsigned long)geo.nDetecU*geo.nDetecV*PROJ_PER_BLOCK);
+                    }
+                }
+            }
+
+            
+            // Now, lets get out the projections from the previous execution of the kernels.
+            if (i>0){
+                for (dev = 0; dev < deviceCount; dev++){
+                    // copy result to host
+                    cudaSetDevice(dev);
+                    cudaMemcpyAsync(result[(i-1)*PROJ_PER_BLOCK+dev*nangles_device], dProjection[(int)(!(i%2))+dev*deviceCount], num_bytes_proj, cudaMemcpyDeviceToHost,stream[dev*2+1]);
+                }
+            }
+            // Make sure Computation on kernels has finished before we launch the next batch.
+            for (dev = 0; dev < deviceCount; dev++){
+                cudaSetDevice(dev);
+                cudaStreamSynchronize(stream[dev*2]);
+            }   
+        }// END angles loop
+        
+        
+        // We still have the last one to get out, do that one
+        int angles_last_device=(nangles-(deviceCount-1)*nangles_device);
+        int size_last_block=nangles_device-(i-1)*PROJ_PER_BLOCK;
+        
+        for (dev = 0; dev < deviceCount; dev++){
+            if(dev+1==deviceCount){
+                size_last_block=angles_last_device-(i-1)*PROJ_PER_BLOCK;
+            }
+            // copy result to host
+            cudaSetDevice(dev);
+            cudaDeviceSynchronize();
+            if(cudaCheckErrors("Fail memcopy fail")){return 1;}
+            cudaMemcpyAsync(result[(i-1)*PROJ_PER_BLOCK+dev*nangles_device], dProjection[(int)(((i-1)%2))+dev*deviceCount], size_last_block*geo.nDetecV*geo.nDetecU*sizeof(float), cudaMemcpyDeviceToHost,stream[dev*2+1]);
+        }
+        // Free memory for the next piece of image
+       
+        cudaDeviceSynchronize();
+    }//END image splits loop
+     
+    return 0;
+    if(cudaCheckErrors("Main loop  fail")){return 1;}
+    ///////////////////////////////////////////////////////////////////////
+    ///////////////////////////////////////////////////////////////////////
+    for (dev = 0; dev < deviceCount; dev++){
+            cudaSetDevice(dev);
+            cudaDestroyTextureObject(texImg[dev]);
+            cudaFreeArray(d_cuArrTex[dev]);
+    }
+    // Freeing Stage
+    for (dev = 0; dev < deviceCount; dev++){
+        cudaSetDevice(dev);
+        cudaFree(dProjection[dev*deviceCount]);
+        cudaFree(dProjection[dev*deviceCount+1]);
+        
+    }
+    free(dProjection);
+    
+    if(!fits_in_memory){
+        for (dev = 0; dev < deviceCount; dev++){
+            cudaSetDevice(dev);
+            cudaFree(dProjection_accum[dev*deviceCount]);
+            cudaFree(dProjection_accum[dev*deviceCount+1]);
+            
+        }
+        free(dProjection_accum);
+    }
+    freeGeoArray(splits,geoArray);
+    cudaFreeHost(projParamsArrayHost);
+   
+    
+    for (int i = 0; i < nStreams; ++i)
+        cudaStreamDestroy(stream[i]) ;
+    
+        if (isHostRegisterSupported & splits>1){
+            cudaHostUnregister(img);
+    }
+    if(cudaCheckErrors("cudaFree  fail")){return 1;}
+    
+//     cudaDeviceReset();
+    return 0;
+}
+
+
+
+
+void CreateTexture(int num_devices,const float* imagedata,Geometry geo,cudaArray** d_cuArrTex, cudaTextureObject_t *texImage,bool alloc)
+{
+    //size_t size_image=geo.nVoxelX*geo.nVoxelY*geo.nVoxelZ;
     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;
-}
-
-
+    if(alloc){
+    for (unsigned int i = 0; i < num_devices; i++){
+        cudaSetDevice(i);
+        
+        //cudaArray Descriptor
+        cudaChannelFormatDesc channelDesc = cudaCreateChannelDesc<float>();
+        //cuda Array
+        cudaMalloc3DArray(&d_cuArrTex[i], &channelDesc, extent);
+    }
+    }
+    for (unsigned int i = 0; i < num_devices; i++){
+        cudaSetDevice(i);
+        cudaMemcpy3DParms copyParams = {0};
+        //Array creation
+        copyParams.srcPtr   = make_cudaPitchedPtr((void *)imagedata, extent.width*sizeof(float), extent.width, extent.height);
+        copyParams.dstArray = d_cuArrTex[i];
+        copyParams.extent   = extent;
+        copyParams.kind     = cudaMemcpyHostToDevice;
+        cudaMemcpy3DAsync(&copyParams);
+    }
+    for (unsigned int i = 0; i < num_devices; i++){
+        cudaSetDevice(i);
+        cudaResourceDesc    texRes;
+        memset(&texRes, 0, sizeof(cudaResourceDesc));
+        texRes.resType = cudaResourceTypeArray;
+        texRes.res.array.array  = d_cuArrTex[i];
+        cudaTextureDesc     texDescr;
+        memset(&texDescr, 0, sizeof(cudaTextureDesc));
+        texDescr.normalizedCoords = false;
+        texDescr.filterMode = cudaFilterModePoint;
+        texDescr.addressMode[0] = cudaAddressModeBorder;
+        texDescr.addressMode[1] = cudaAddressModeBorder;
+        texDescr.addressMode[2] = cudaAddressModeBorder;
+        texDescr.readMode = cudaReadModeElementType;
+        cudaCreateTextureObject(&texImage[i], &texRes, &texDescr, NULL);
+        
+    }
+    for (unsigned int i = 0; i < num_devices; i++){
+        cudaSetDevice(i);
+        cudaDeviceSynchronize();
+    }
+    cudaCheckErrors("Texture object creation fail");
+}
+
+/* This code generates the geometries needed to split the image properly in
+ * cases where the entire image does not fit in the memory of the GPU
+ **/
+void splitImage(unsigned int splits,Geometry geo,Geometry* geoArray, unsigned int nangles){
+    
+    unsigned long splitsize=(geo.nVoxelZ+splits-1)/splits;// ceil if not divisible
+    for(unsigned int sp=0;sp<splits;sp++){
+        geoArray[sp]=geo;
+        // All of them are splitsize, but the last one, possible
+        geoArray[sp].nVoxelZ=((sp+1)*splitsize<geo.nVoxelZ)?  splitsize:  geo.nVoxelZ-splitsize*sp;
+        geoArray[sp].sVoxelZ= geoArray[sp].nVoxelZ* geoArray[sp].dVoxelZ;
+        
+        // We need to redefine the offsets, as now each subimage is not aligned in the origin.
+        geoArray[sp].offOrigZ=(float *)malloc(nangles*sizeof(float));
+        for (unsigned int i=0;i<nangles;i++){
+            geoArray[sp].offOrigZ[i]=geo.offOrigZ[i]-geo.sVoxelZ/2+sp*geoArray[0].sVoxelZ+geoArray[sp].sVoxelZ/2;
+        }
+        
+    }
+    
+}
 
 /* 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 <math.h>
#include <stdio.h>
#include "projection.hpp"

inline int cudaCheckErrors(const char * msg)
{
   cudaError_t __err = cudaGetLastError();
   if (__err != cudaSuccess)
   {
      printf("CUDA:Siddon_projection:%s:%s\n",msg, cudaGetErrorString(__err));
      cudaDeviceReset();
      return 1;
   }
   return 0;
}
    
    
// Declare the texture reference.
    texture<float, cudaTextureType3D , cudaReadModeElementType> tex;
#define MAXTREADS 1024
/*GEOMETRY DEFINITION
 *
 *                Detector plane, behind
 *            |-----------------------------|
 *            |                             |
 *            |                             |
 *            |                             |
 *            |                             |
 *            |      +--------+             |
 *            |     /        /|             |
 *   A Z      |    /        / |*D           |
 *   |        |   +--------+  |             |
 *   |        |   |        |  |             |
 *   |        |   |     *O |  +             |
 *    --->y   |   |        | /              |
 *  /         |   |        |/               |
 * V X        |   +--------+                |
 *            |-----------------------------|
 *
 *           *S
 *
 *
 *
 *
 *
 **/



__global__ void kernelPixelDetector( Geometry geo,
        float* detector,
        Point3D source ,
        Point3D deltaU,
        Point3D deltaV,
        Point3D uvOrigin){
    
//     size_t idx = threadIdx.x + blockIdx.x * blockDim.x;
    
    unsigned long y = blockIdx.y * blockDim.y + threadIdx.y;
    unsigned long x = blockIdx.x * blockDim.x + threadIdx.x;
    unsigned long idx =  x  * geo.nDetecV + y;
    
    if ((x>= geo.nDetecU) | (y>= geo.nDetecV))
        return;
    
    
    
    
    /////// Get coordinates XYZ of pixel UV
    int pixelV = geo.nDetecV-y-1;

    float axm,aym,azm;
    float axM,ayM,azM;
    // In the paper Nx= number of X planes-> Nvoxel+1
   
    axm=fminf(__fdividef(-source.x,ray.x),__fdividef(geo.nVoxelX-source.x,ray.x));
    aym=fminf(__fdividef(-source.y,ray.y),__fdividef(geo.nVoxelY-source.y,ray.y));
    azm=fminf(__fdividef(-source.z,ray.z),__fdividef(geo.nVoxelZ-source.z,ray.z));
    axM=fmaxf(__fdividef(-source.x,ray.x),__fdividef(geo.nVoxelX-source.x,ray.x));
    ayM=fmaxf(__fdividef(-source.y,ray.y),__fdividef(geo.nVoxelY-source.y,ray.y));
    azM=fmaxf(__fdividef(-source.z,ray.z),__fdividef(geo.nVoxelZ-source.z,ray.z));
    
    float am=fmaxf(fmaxf(axm,aym),azm);
    float aM=fminf(fminf(axM,ayM),azM);
    
    // line intersects voxel space ->   am<aM
    if (am>=aM)
