# Macro to compute the RMSE between the Shepp-Logan head phantom
# and a processed version (e.g., after tomography reconstruction)

macro rmse:

	gs = generate ("head", 512, 1);		# Gold-standard, orig. headphantom
	cmp = mainimage;					# Compare to main image
	a = cmp-gs;
	a = a*a;							# Squared deviation
	mse = measure (a, 0);				# Mean value of a is MSE
	print ("Root mean squared error is ", power(mse, 0.5));

endmacro;

#################################################################

procedure numproj:					# !! disable this line and...
#macro numproj:						# !! enable this line to enable the macro

	hp = generate ("head", 512, 1);		# Create a head phantom
	dim (anginc[5], result[5]);
	anginc[0]=20; anginc[1]=10; anginc[2]=5;
	anginc[3]=2; anginc[4]=1;

	for (i=0, 4);

		sino = sinogram (hp, int(anginc[i]));		# Create sinogram
		#display (sino);
		rec = reconstruct (sino, "fbp360", 0,0.5,0);		# Rewconstruction w/FBP, default options
		display (rec);

		a = hp - rec;
		a = a*a;
		mse = measure (a, 0);				# Mean value of a is MSE
		j = 360/anginc[i];
		print ("For ", j, " views, RMSE is ", power(mse, 0.5));
		anginc[i] = j;
		result[i] = power(mse, 0.5);
	endfor;
	graph (4, anginc, result);

endproc;			# !! Also adjust this to "endmacro"