About optional parameters <init = "multigrid">

[Kwong-Yam-Lie]

Hi, What does this initialization(init = "multigrid") do? I tried the case with this parameter, but it's too slow...

imgsart = algs.ossart(prosart,x.geo,x.angles,niter = 1, init = "multigrid", verbose = False)

I have read related source code.

TIGRE/Python/tigre/utilities/init_multigrid.py

Lines 7 to 52 in c55a6bb

	def init_multigrid(proj, geo, alpha, alg):
	# WARNING: This takes a lot of memory!
	if alg == "SART":
	italg = tigre.algorithms.sart
	if alg == "SIRT":
	italg = tigre.algorithms.sirt
	finalsize = geo.nVoxel
	maxval = max(proj.ravel())
	minval = min(proj.ravel())
	# Start with 16 (raise this for larger images)
	geo.nVoxel = np.array([16, 16, 16])
	geo.dVoxel = geo.sVoxel / geo.nVoxel
	if (geo.nVoxel > finalsize).all():
	return np.zeros(finalsize, dtype=np.float32)
	niter = 100
	initres = np.zeros(geo.nVoxel, dtype=np.float32)
	while (geo.nVoxel != finalsize).all():
	geo.dVoxel = geo.sVoxel / geo.nVoxel
	initres = italg(proj, geo, alpha, niter, init=initres, verbose=False)
	# get new dims(should be a way to do this more efficiently).
	geo.nVoxel = geo.nVoxel * 2
	geo.nVoxel[geo.nVoxel > finalsize] = finalsize[geo.nVoxel > finalsize]
	geo.dVoxel = geo.sVoxel / geo.nVoxel
	(x, y, z) = (
	np.linspace(minval, maxval, geo.nVoxel[0] / 2, dtype=np.float32),
	np.linspace(minval, maxval, geo.nVoxel[1] / 2, dtype=np.float32),
	np.linspace(minval, maxval, geo.nVoxel[2] / 2, dtype=np.float32),
	)
	# evaluate the function sart at the points xv,yv,zv
	xv, yv, zv = [
	tile_array(tile_array(x, 2), geo.nVoxel[0] ** 2),
	tile_array(tile_array(y, 2), geo.nVoxel[0] ** 2),
	tile_array(tile_array(x, 2), geo.nVoxel[0] ** 2),
	]
	initres = RegularGridInterpolator((x, y, z), initres)(np.column_stack((xv, yv, zv)))
	initres = initres.reshape(geo.nVoxel)
	return initres

Is the number of iterations set to 100 (niter = 100)too large?

TIGRE/Python/tigre/utilities/init_multigrid.py

Line 23 in c55a6bb

niter = 100

I'm curious what's the difference between the <init = "multigrid"> and the <init = "FDK">.
In theory, do these two iteration methods have any effect on iteration speed and accuracy?

Can you give me a brief explanation or recommend related articles？ Thanks！

[AnderBiguri]

Hi @Kwong-Yam-Lie
Admitedly I have not used multigrid a lot, it was an idea from quite long ago that has not been explored much.
The idea is to recosntruct in a smaller image size (bigger voxel size) in steps until you have the size of the real problem. I agree, 100 seems too large for the initialization.

FDK initializes the algorithm with the FDK image as initial, instead of all zeroes. This is useful when you don't have a ton of time for the recon, but often the FDK artifacts are still there after several iterations, if initialized like this.

[Kwong-Yam-Lie]

I understand. At first, I thought you were doing this to speed up the reconstruction or improve the accuracy of the reconstruction results, maybe there is some mathematical conclusion that I don't know. But after a simple verification, I found that there was no good result.

I focused on this mainly because I had an idea: we could reconstruct pixels in the center of an object while maintaining voxel size, and then gradually spread towards the object boundary. Equivalently, we first take the local data of the detector and then extend it to the global data of the probe.

This is similar to what you think, but it's not clear to me if you've already done this in backprojection (Atb). If you have already tried, can you tell me the feasibility of this scheme?

[AnderBiguri]

@Kwong-Yam-Lie In theory, it should speed up recon, in the sense that you'd need less iterations to get to the same results, because your initial image should be a good approximation to the solution.

Your idea, if I understand correctly, is to have some sort of mixed-scale mesh to recosntruct, where pixels are small in the center and bigger outside? Or something else?

Atb is quite simple, it just updates each voxel with the information in the detector in the corresponding location. The pixels used are the same as in the desired image, nothing else.

[Kwong-Yam-Lie]

In my idea, pixels are equal in center and outside.
It's just that each iteration starts with some rays from the center of the detector, and then uses more until all rays are used. The purpose is to try to speed up the convergence of pixels outside.

My problem has been solved, thanks!