             ac=az;
             az+=azu;
         }
-        aminc=min(min(ax,ay),az);
+        aminc=fminf(fminf(ax,ay),az);
     }
     detector[idx]=sum*maxlength;
 }

        
        for (unsigned int i=0;i<nangles;i+=(unsigned int)deviceCount){
            for (dev = 0; dev < deviceCount; dev++){
                geoArray[sp].alpha=angles[(i+dev)*3];
                geoArray[sp].theta=angles[(i+dev)*3+1];
                geoArray[sp].psi  =angles[(i+dev)*3+2];
                //precomute distances for faster execution
                //Precompute per angle constant stuff for speed
                computeDeltas_Siddon(geoArray[sp],i+dev, &uvOrigin, &deltaU, &deltaV, &source);
                //Ray tracing!
               
                cudaSetDevice(dev);
                kernelPixelDetector<<<grid,block>>>(geoArray[sp],dProjection[dev], source, deltaU, deltaV, uvOrigin,texImg[dev]);
                //cudaCheckErrors("Kernel fail");
                
            }
            
            for (dev = 0; dev < deviceCount; dev++){
                cudaSetDevice(dev);
                if (!fits_in_memory){
                    mexPrintf("bipbop %d\n",(int)sp);
                    cudaMemcpyAsync(dProjection_accum[dev], result[i+dev], num_bytes_proj, cudaMemcpyHostToDevice);
                    vecAddInPlace<<<(geo.nDetecU*geo.nDetecV+MAXTREADS-1)/MAXTREADS,MAXTREADS>>>(dProjection[dev],dProjection_accum[dev],(unsigned long)geo.nDetecU*geo.nDetecV );
                }
            }
            for (dev = 0; dev < deviceCount; dev++){
                // copy result to host
                cudaSetDevice(dev);
                cudaMemcpyAsync(result[i+dev], dProjection[dev], num_bytes_proj, cudaMemcpyDeviceToHost);
                
            }
            cudaCheckErrors("cudaMemcpy output fail");

    float imin,imax,jmin,jmax,kmin,kmax;
    // for X
    if( source.x<pixel1D.x){
        imin=(am==axm)? 1             : ceilf (source.x+am*ray.x);
        imax=(aM==axM)? geo.nVoxelX   : floorf(source.x+aM*ray.x);
    }else{
        imax=(am==axm)? geo.nVoxelX-1 : floorf(source.x+am*ray.x);
        imin=(aM==axM)? 0             : ceilf (source.x+aM*ray.x);
    }
    // for Y
    if( source.y<pixel1D.y){
        jmin=(am==aym)? 1             : ceilf (source.y+am*ray.y);
        jmax=(aM==ayM)? geo.nVoxelY   : floorf(source.y+aM*ray.y);
    }else{
        jmax=(am==aym)? geo.nVoxelY-1 : floorf(source.y+am*ray.y);
        jmin=(aM==ayM)? 0             : ceilf (source.y+aM*ray.y);
    }
    // for Z
    if( source.z<pixel1D.z){
        kmin=(am==azm)? 1             : ceilf (source.z+am*ray.z);
        kmax=(aM==azM)? geo.nVoxelZ   : floorf(source.z+aM*ray.z);
    }else{
        kmax=(am==azm)? geo.nVoxelZ-1 : floorf(source.z+am*ray.z);
        kmin=(aM==azM)? 0             : ceilf (source.z+aM*ray.z);
    }
    
    // get intersection point N1. eq(20-21) [(also eq 9-10)]
