What is (or what are) the best conventional method(s)/algorithm(s) to form clusters of active t-value assigned voxels from fMRI results in Matlab?

I have imported into Matlab an excel sheet with over 500,000 voxels each assigned a t value and x,y,z coordinates based on fMRI BOLD signal results. I would like to form clusters of activation in 6 regions of the brain from the activated voxels. I would like to know the best way (possibly method or algorithm) to form clusters (accounting for overlapping active voxels) in Matlab, so that I am able to calculate the number of active clusters and volumetric features, such as volume of the clusters. I know one method is a density-based method which could scan a cube of 3*3*3 region and consider that region to be a cluster if a certain number of active voxels exist in there, but I want to know which method(s) in matlab is best to form clusters so that they best include active voxels regardless of shape.