Commits

kritisen committed 30e2c2a

FBP short scan capability added as M file. SART modified to be updated only over circular object support, and cleaned up FBP code.

- FBP short scan adjustments done by MATLAB function: weightForShortScan()
- FBP code cleaned up
- SART: added capability to consider only circular object support
- SART: adjustments within circular object support weighted by Hamming window as advised in Kak book

Comments (0)

Files changed (5)

 
 	return intsN;
 }
+
+//same as Prj_getPrjCoefByDetr(..) except that only pixels lying within object support are considered
+// size of object support is defined by radius r0 from center of reconstruction area
+int Prj_getPrjCoefByDetr_ObjSupport(vertex p, vector n, double *xs, double dx, int nxs, double *ys, 
+			double dy, int nys, double ddet, double *coef, int *offset, double * sumc, double r0)
+{
+	double kx = 1/dx;
+	double ky = 1/dy;
+
+	vertex s0,s1;
+	vertex s, snextc, snextr;
+	vector ns;
+	double tc,tr,c,sumcoef=0;
+	double xl,xr,yd,yu;
+
+	int dr, dc, c_off, r_off, i_off;
+
+	int intsN=0;
+	double prj = 0;
+	double w;
+
+	//scaling
+	p.x = kx * p.x;
+	p.y = ky * p.y;
+	n.x = kx * n.x;
+	n.y = ky * n.y;
+	
+    w = 1;
+    xl = *xs;
+	xr = *(xs+nxs-1);
+	yd = *ys;
+	yu = *(ys+nys-1);
+	
+	if(abs(n.x) >= abs(n.y)){
+		
+		if(Prj_get_s0s1_x(p.x, p.y, n.x, n.y, xl, xr, yd, yu, &s0, &s1, &ns)!=0)
+			return -1;
+
+		if(s0.y<=s1.y){
+
+			s=s0;
+
+			c_off = floor(s.x-*xs);
+			snextc.x = floor(s.x+1);
+			tc = (snextc.x-s.x)/ns.x;
+			snextc.y = ns.y * tc + s.y;
+
+			r_off = floor(s.y-*ys);
+			snextr.y = floor(s.y+1);
+			tr = (snextr.y-s.y)/ns.y;
+			snextr.x = ns.x * tr + s.x;
+
+			while(s.x<s1.x-1e-8){
+				if(tc<tr){					
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tc * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+
+					//update s, snextc
+					s = snextc;
+					tr = tr - tc;
+					c_off ++;
+					snextc.x = s.x+1;
+					tc = 1/ns.x;
+					snextc.y = ns.y * tc + s.y;
+				}
+				else if(tc>tr){
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tr * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+
+					//update s, snextr
+					s = snextr;
+					tc = tc - tr;
+					r_off ++;
+					snextr.y = s.y+1;
+					tr = 1/ns.y;
+					snextr.x = ns.x * tr + s.x;
+				}	
+				else{ //snextc.x == snextr.x
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tc * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+
+					//update s, snextr, snextc
+					s = snextc;
+					r_off ++;
+					snextr.y = s.y+1;
+					tr = 1/ns.y;
+					snextr.x = ns.x * tr + s.x;
+					c_off ++;
+					snextc.x = s.x+1;
+					tc = 1/ns.x;
+					snextc.y = ns.y * tc + s.y;
+				}
+			}
+
+
+		}
+		else{ //s0.y > s1.y
+
+			s=s0;
+
+			c_off = floor(s.x-*xs);
+			snextc.x = floor(s.x+1);
+			tc = (snextc.x-s.x)/ns.x;
+			snextc.y = ns.y * tc + s.y;
+
+			r_off = ceil(s.y-1-*ys);
+			snextr.y = ceil(s.y-1);
+			tr = (snextr.y-s.y)/ns.y;
+			snextr.x = ns.x * tr + s.x;
+
+			while(s.x<s1.x-1e-8){
+				if(tc<tr){					
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tc * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+					
+					//update s, snextc
+					s = snextc;
+					tr = tr - tc;
+					c_off ++;
+					snextc.x = s.x+1;
+					tc = 1/ns.x;
+					snextc.y = ns.y * tc + s.y;
+				}
+				else if(tc>tr){
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tr * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+
+					//update s, snextr
+					s = snextr;
+					tc = tc - tr;
+					r_off --;
+					snextr.y = s.y-1;
+					tr = -1/ns.y;
+					snextr.x = ns.x * tr + s.x;
+				}	
+				else{ //snextc.x == snextr.x
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tc * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+
+					//update s, snextr, snextc
+					s = snextc;
+					c_off ++;
+					snextc.x = s.x+1;
+					tc = 1/ns.x;
+					snextc.y = ns.y * tc + s.y;
+					r_off --;
+					snextr.y = s.y-1;
+					tr = -1/ns.y;
+					snextr.x = ns.x * tr + s.x;
+				}
+			}
+		}
+		
+	}
+	else{
+
+		if(Prj_get_s0s1_y(p.x, p.y, n.x, n.y, xl, xr, yd, yu, &s0, &s1, &ns)!=0)
+			return -1;
+		
+		if(s0.x<=s1.x){
+
+			s=s0;
+
+			r_off = floor(s.y-*ys);
+			snextr.y = floor(s.y+1);
+			tr = (snextr.y-s.y)/ns.y;
+			snextr.x = ns.x * tr + s.x;
+
+			c_off = floor(s.x-*xs);
+			snextc.x = floor(s.x+1);
+			tc = (snextc.x-s.x)/ns.x;
+			snextc.y = ns.y * tc + s.y;
+
+			while(s.y<s1.y-1e-8){
+				if(tr<tc){					
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tr * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+					
+					//update s, snextr
+					s = snextr;
+					tc = tc - tr;
+					r_off ++;
+					snextr.y = s.y+1;
+					tr = 1/ns.y;
+					snextr.x = ns.x * tr + s.x;
+				}
+				else if(tr>tc){
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tc * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+
+					//update s, snextr
+					s = snextc;
+					tr = tr - tc;
+					c_off ++;
+					snextc.x = s.x+1;
+					tc = 1/ns.x;
+					snextc.y = ns.y * tc + s.y;
+				}	
+				else{ //snextr.y == snextc.y
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tr * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+
+					//update s, snextr, snextc
+					s = snextr;
+					c_off ++;
+					snextc.x = s.x+1;
+					tc = 1/ns.x;
+					snextc.y = ns.y * tc + s.y;
+					r_off ++;
+					snextr.y = s.y+1;
+					tr = 1/ns.y;
+					snextr.x = ns.x * tr + s.x;
+				}
+			}
+		}
+		else{ //s0.x > s1.x
+
+			s=s0;
+
+			r_off = floor(s.y-*ys);
+			snextr.y = floor(s.y+1);
+			tr = (snextr.y-s.y)/ns.y;
+			snextr.x = ns.x * tr + s.x;
+
+			c_off = ceil(s.x-1-*xs);
+			snextc.x = ceil(s.x-1);
+			tc = (snextc.x-s.x)/ns.x;
+			snextc.y = ns.y * tc + s.y;
+
+			while(s.y<s1.y-1e-8){
+				if(tr<tc){					
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tr * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+					
+					//update s, snextc
+					s = snextr;
+					tc = tc - tr;
+					r_off ++;
+					snextr.y = s.y+1;
+					tr = 1/ns.y;
+					snextr.x = ns.x * tr + s.x;
+				}
+				else if(tr>tc){
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tc * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+
+					//update s, snextr
+					s = snextc;
+					tr = tr - tc;
+					c_off --;
+					snextc.x = s.x-1;
+					tc = -1/ns.x;
+					snextc.y = ns.y * tc + s.y;
+				}	
+				else{ //snextc.x == snextr.x
+					//save
+					double dist = sqrt(pow((nys-2-r_off) - (nys-1)/2, 2) + pow(c_off - (nxs-1)/2, 2));
+					if (dist <= r0) {	//if pixel lies within object support
+						c = tr * w;
+						*coef++ = c;
+						sumcoef += c;
+						*offset++ = (nys-1)*c_off + (nys-2-r_off);
+						intsN++;
+					}
+
+					//update s, snextr, snextc
+					s = snextr;
+					r_off ++;
+					snextr.y = s.y+1;
+					tr = 1/ns.y;
+					snextr.x = ns.x * tr + s.x;
+					c_off --;
+					snextc.x = s.x-1;
+					tc = -1/ns.x;
+					snextc.y = ns.y * tc + s.y;
+				}
+			}
+		}
+	}
+
+	if(offset)
+		*offset = -1;
+
+	if(sumc)
+		*sumc = sumcoef;
+
+	return intsN;
+}
 
 extern int Prj_getInts(vertex p, vector n, double *xs, double dx, int nxs, double *ys, double dy, int nys, vertex *ints);
 extern int Prj_getPrjCoefByDetr(vertex p, vector n, double *xs, double dx, int nxs, double *ys, double dy, int nys, double ddet, double *coef, int *offset, double * sumc);
+extern int Prj_getPrjCoefByDetr_ObjSupport(vertex p, vector n, double *xs, double dx, int nxs, double *ys, double dy, int nys, double ddet, double *coef, int *offset, double * sumc, double r0);
 extern double Prj_getPrjByDetr(vertex p, vector n, double *img, double *xs, double dx, int nxs, double *ys, double dy, int nys, double ddet, double *coef, int *offset, double * sumc);
 
 #endif	// Prj_h_
     if(dsamp > 1)
         mexPrintf("Use %dX Downsampling.\n",dsamp);mexEvalString("drawnow;");
     
-    //allocate memory space needed for single projection ray
-    //double * pcoef = (double *)calloc(imgN,sizeof(double));
-    //int * poffset = (int *)calloc(imgN,sizeof(int));
-    
 	mexPrintf("\nCenter shift =%lf\n", ra.cShift);mexEvalString("drawnow;");
     for(int angi=0; angi<pa.angN; angi++){
 		if(angi%10==0){ mexPrintf(".");mexEvalString("drawnow;"); }
 			if (ia.pxlW == 1 && ddet == 1){
 				double D = pa.scanR;	//distance from source to object center 
 				double tau = ra.cShift;	//center shift
+				
+				//double gammaM = atan((pa.detN/2)/D); //weighting to adjust for short-scan
+				//for(int deti=0; deti<pa.detN; deti++){
+				//	double s = (deti - pa.detN/2);  
+				//	double gamma = atan(s/D);
+				//	double L2 = 2*(gammaM - gamma);
+				//	double L3 = pi - 2 * gamma;
+				//	double L4 = pi + 2 * gammaM;
+				//	double Ea = pow(sin(pi/4*ang/(gammaM - gamma)), 2);
+				//	double Ec = pow(sin(pi/4*(pi+2*gammaM - ang)/(gammaM + gamma)), 2);
+				//	double factor = 0;	// weight the projections which lie outside this range 
+										// 	with the factor 0. so that they do not have any effect
+										//	in the FBP
+				//	if (ang >= 0 && ang <= L2)
+				//		factor = Ea;
+				//	else if (ang >= L2 && ang <= L3)
+				//		factor = 1;
+				//	else if (ang >= L3 && ang <= L4)
+				//		factor = Ec;
+				//	prj[pa.detN*angi+deti] *= factor;
+				//} // **** This is carried out in MATLAB function >>weightForShortScan()<< cefore calling fanrec() 
+				//			This accounts for short scans which do not cover the whole angular range 2pi
+				//				but run through pi+2gammaM	***
+				
 				//for(int deti=0; deti<pa.detN; deti++){
 				//	double s = (deti - pa.detN/2);     
 				//	prj[pa.detN*angi+deti] *= (D-(tau*s)/D)/(sqrt(D*D+s*s));
-				//} // **** This is carried out in MATLAB cefore calling fanrec() ***
+				//} // **** This is carried out in MATLAB function >>weightFanProjection()<< 
+				//				before calling fanrec() 
+				//			This function accounts for different FBP formula for fanbeam than parallel beam
+				//				and also accounts for center shift***
 				 
 				//filterProjection(prj + pa.detN*angi, pa.detN);	//filter the projection with the ramp 
 																//filter h(s) / smoothed ramp filter
 				// **** filterProjection() is carried out in MATLAB cefore calling fanrec(),
 				// by calling MATLAB code for filtering projections ***
 				
-				//if (floor(ang*180/pi) == 0){ //for debug
-					//mexPrintf("\nang =%lf\n", ang);mexEvalString("drawnow;");
 				for (int j = 0; j < imgN; j++) //loop through all image pixels
 				{
 					//calculate {r, phi} for (x, y) with origin as (0, tau)
 						}
 					} //if sDash lies within data range
 				} //end loop through all image pixels
-				//}
-			} 
+			} //check that pixel size is 1 and detector-RA distance is zero
 			else {
 				mexPrintf("FBP error: pixel size other than 1, and non-zero detector distance is not yet supported. Best to rescale all distances to image pixel size, and use detector to rotation axis distance as zero.");mexEvalString("drawnow;");
 				return -1;		
     for(int i=0;i<imgN;i++)
         img[i]*=da;
     
-    ////do Filtration
-    //if(hRec->flt && hImg->fa.type)
-	//	Img_doFlt(hImg);
-
-	//free(pcoef);
-	//free(poffset);
-    
     mexPrintf("Done\n");mexEvalString("drawnow;");
 
     return 0;      
     double * pcoef = (double *)calloc(imgN,sizeof(double));
     int * poffset = (int *)calloc(imgN,sizeof(int));
     	
-	mexPrintf("\nCenter shift =%lf\n", ra.cShift);mexEvalString("drawnow;");
+	if (ra.cShift != 0) { mexPrintf("\nCenter shift =%lf\n", ra.cShift);mexEvalString("drawnow;"); }
 	for(int iteri=0;iteri<ra.iterN;iteri++)//iteration loop
 	{
         mexPrintf("\niter: %d ",iteri);mexEvalString("drawnow;");
 						//calculate new Proj
 						double sumc=0, dprj=0;
 		
-						int n = Prj_getPrjCoefByDetr(pRay, nRay, ia.xs, ia.pxlW, ia.colN+1, ia.ys, ia.pxlH, ia.rowN+1, ddet, pcoef, poffset, &sumc);
+						double r0 = pa.detN/2 * pa.scanR / sqrt(pow(pa.scanR, 2) + pow(pa.detN/2, 2));
+						int n = Prj_getPrjCoefByDetr_ObjSupport(pRay, nRay, ia.xs, ia.pxlW, ia.colN+1, ia.ys, 
+														ia.pxlH, ia.rowN+1, ddet, pcoef, poffset, &sumc, r0);
 						if(n){
 							for(int i=0; i<n; i++){
 								dprj += img[poffset[i]] * pcoef[i];
 						if(dprj>=0){
 							dprj = prj[angidx*pa.detN+deti] - dprj;
 							for(int i=0;i<n;i++){
-								dimg[poffset[i]] += pcoef[i] * dprj / sumc / dsamp;
+								double hammWeight = 0.5 - 0.5 * cos(2*pi*i/n);	//0.54 - 0.46 * cos(2*pi*i/n);
+								dimg[poffset[i]] += hammWeight * pcoef[i] * dprj / sumc / dsamp;
 								sumDet[poffset[i]] += pcoef[i] / dsamp;
 							}
 						}

utilities/weightForShortScan.m

+function [ P2 ] = weightForShortScan( P, D, theta )
+%weightForShortScan - weight the projections for short scan
+%   according to Kak and Slaney Section 3.5, pg 99
+%   P - sinogram
+%   D - source to rotation axis distance (scaled by pixel size)
+%   theta - angles corresponding to sinogram 
+
+[detN nAng] = size(P);
+P2 = P;
+gammaM = atan((detN/2)/D);
+for angi = 1: nAng
+    ang = theta(angi);
+    ang2 = ang - theta(1); % formula based on angles between [0, pi + 2*gammaM]
+                            % but Xradia scanner acquires:
+                            % [-pi/2-gammaM, pi/2 + gammaM]
+    for deti = 1: detN
+        s = (deti-1) - detN/2;
+        gamma = atan(s/D);
+        L2 = 2*(gammaM - gamma);
+        L3 = pi - 2 * gamma;
+        L4 = pi + 2 * gammaM;
+        Ea = sin(pi/4*ang2/(gammaM - gamma))^2;
+        Ec = sin(pi/4*(pi+2*gammaM - ang2)/(gammaM + gamma))^2;
+        factor = 0;	% weight the projections which lie outside this range 
+                        % 	with the factor 0. so that they do not have any effect
+                        %	in the FBP
+        if (ang2 >= 0 && ang2 <= L2)
+            factor = Ea;
+        elseif (ang2 >= L2 && ang2 <= L3)
+            factor = 1;
+        elseif (ang2 >= L3 && ang2 <= L4)
+            factor = Ec;
+        end
+        P2(deti, angi) = P2(deti, angi) * factor;
+    end
+end
+
+end
+