From 39904689f4dfbea5f5e00917959f720244f407bd Mon Sep 17 00:00:00 2001
From: AnderBiguri <aabhca20@EEPC-301.campus.bath.ac.uk>
Date: Tue, 7 Feb 2017 17:58:45 +0000
Subject: [PATCH] Fix bug in voxel_backprojection2

CGLS was behaving weirdly as there was a weight missing.
---
 Algorithms/CGLS.m               |  4 +++-
 Demos/d09_Algorithms04.m        |  2 +-
 Source/voxel_backprojection2.cu | 18 ++++++++++++++----
 3 files changed, 18 insertions(+), 6 deletions(-)

diff --git a/Algorithms/CGLS.m b/Algorithms/CGLS.m
index 0e7e3a98..f1f5742b 100644
--- a/Algorithms/CGLS.m
+++ b/Algorithms/CGLS.m
@@ -141,7 +141,9 @@
     if ii>1 && errorL2(ii)>errorL2(ii-1)
         % OUT!
        x=x-alpha*p;
-
+       if verbose
+          disp(['CGLS stoped in iteration N', num2str(ii),' due to divergence.']) 
+       end
        return; 
     end
     % If step is adecuatem, then continue withg CGLS
diff --git a/Demos/d09_Algorithms04.m b/Demos/d09_Algorithms04.m
index 10b4a978..4762f685 100644
--- a/Demos/d09_Algorithms04.m
+++ b/Demos/d09_Algorithms04.m
@@ -108,7 +108,7 @@
 
 
 %   'alpha_red':   Defines the reduction rate of the TV hyperparameter
-alpha_red=0.95;
+alpha_red=0.noise_projections95;
 
 %   'Ratio':       The maximum allowed image/TV update ration. If the TV 
 %                  update changes the image more than this, the parameter
diff --git a/Source/voxel_backprojection2.cu b/Source/voxel_backprojection2.cu
index b5f82b7d..664d4347 100644
--- a/Source/voxel_backprojection2.cu
+++ b/Source/voxel_backprojection2.cu
@@ -90,7 +90,13 @@ do { \
      *
      **/
     texture<float, cudaTextureType3D , cudaReadModeElementType> tex;
-
+__global__ void matrixConstantMultiply(const Geometry geo,float* image,float constant){
+    size_t idx = threadIdx.x + blockIdx.x * blockDim.x;
+     for(; idx<geo.nVoxelX* geo.nVoxelY *geo.nVoxelZ; idx+=gridDim.x*blockDim.x) {
+            image[idx]*=constant;
+     }
+    
+}
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 // RB, 10/31/2016: Add constant memory arrays to store parameters for all projections to be analyzed during a single kernel call
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
@@ -259,14 +265,16 @@ __global__ void kernelPixelBackprojection(const Geometry geo, float* image,const
             realD.y=-realDaux.x*sinalpha  + realDaux.y*cosalpha; //sin(-x)=-sin(x) , cos(-x)=cos(x)
             float L,l;
             L = sqrt( (realS.x-realD.x)*(realS.x-realD.x)+ (realS.y-realD.y)*(realS.y-realD.y)+ (realD.z)*(realD.z)); // Sz=0 always.
-            l = sqrt( (realS.x-realvoxel.x)*(realS.x-realvoxel.x)+ (realS.y-realvoxel.y)*(realS.y-realvoxel.y)+ (realS.z-realvoxel.z)*(realS.z-realvoxel.z));
+            l = sqrt( (realS.x-realvoxel.x)*(realS.x-realvoxel.x)
+                    + (realS.y-realvoxel.y)*(realS.y-realvoxel.y)
+                    + (realS.z-realvoxel.z)*(realS.z-realvoxel.z));
             weigth=L*L*L/(geo.DSD*l*l);
             
             // Get Value in the computed (U,V) and multiply by the corresponding weigth.
             // indAlpha is the ABSOLUTE number of projection in the projection array (NOT the current number of projection set!)
             voxelColumn[colIdx]+=tex3D(tex, u +0.5 ,
                     v +0.5 ,
-                    indAlpha+0.5)*weigth;
+                    indAlpha+0.5)* weigth;
         }  // END iterating through column of voxels
         
     }  // END iterating through multiple projections
@@ -415,9 +423,11 @@ int voxel_backprojection2(float const * const projections, Geometry geo, float*
         cudaMemcpyToSymbol(projParamsArray2Dev, projParamsArray2Host, sizeof(Point3D)*7*PROJ_PER_KERNEL);
         
         kernelPixelBackprojection<<<grid,block>>>(geo,dimage,i,nalpha);
+
         cudaCheckErrors("Kernel fail");
     }  // END for
-    
+    matrixConstantMultiply<<<60,MAXTREADS>>>( geo,dimage,geo.dVoxelX*geo.dVoxelY*geo.dVoxelZ/(geo.dDetecU*geo.dDetecV));
+
     //////////////////////////////////////////////////////////////////////////////////////
     // END Main reconstruction loop: go through projections (rotation angles) and backproject
     //////////////////////////////////////////////////////////////////////////////////////