     float axm,aym,azm;
     float axM,ayM,azM;
     // In the paper Nx= number of X planes-> Nvoxel+1
-    axm=min(-source.x/ray.x,(geo.nVoxelX-source.x)/ray.x);
-    aym=min(-source.y/ray.y,(geo.nVoxelY-source.y)/ray.y);
-    azm=min(-source.z/ray.z,(geo.nVoxelZ-source.z)/ray.z);
-    axM=max(-source.x/ray.x,(geo.nVoxelX-source.x)/ray.x);
-    ayM=max(-source.y/ray.y,(geo.nVoxelY-source.y)/ray.y);
-    azM=max(-source.z/ray.z,(geo.nVoxelZ-source.z)/ray.z);
-    float am=max(max(axm,aym),azm);
-    float aM=min(min(axM,ayM),azM);
+    axm=fminf(-source.x/ray.x,(geo.nVoxelX-source.x)/ray.x);
+    aym=fminf(-source.y/ray.y,(geo.nVoxelY-source.y)/ray.y);
+    azm=fminf(-source.z/ray.z,(geo.nVoxelZ-source.z)/ray.z);
+    axM=fmaxf(-source.x/ray.x,(geo.nVoxelX-source.x)/ray.x);
+    ayM=fmaxf(-source.y/ray.y,(geo.nVoxelY-source.y)/ray.y);
+    azM=fmaxf(-source.z/ray.z,(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)

        cudaCheckErrors("cudaMalloc projections fail");
    }
    
       

    for (unsigned int sp=0;sp<splits;sp++){
        
        // Create texture objects for all GPUs

    
    // get index of first intersection. eq (26) and (19)
    int i,j,k;
    float aminc=fminf(fminf(ax,ay),az);
    i=(int)floorf(source.x+ (aminc+am)/2*ray.x);
    j=(int)floorf(source.y+ (aminc+am)/2*ray.y);
    k=(int)floorf(source.z+ (aminc+am)/2*ray.z);
    // Initialize
    float ac=am;
    //eq (28), unit anlges
    float axu,ayu,azu;
    axu=1/fabs(ray.x);
    ayu=1/fabs(ray.y);
    azu=1/fabs(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=sqrtf(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
