An incremental form interpolation model integrated with the anisotropic Smolyak method and the Non-Uniform Rational B-Spline (NURBS) method is proposed for multiple material topology optimization. Only a single category of design variables is exploited to characterize the material properties for any number of material types and the Smolyak method embedded within the NURBS method is devised to represent the topology density field while keeping the number of polynomials to a minimum, leading to a significant reduction in the number of design variables. An arbitrary number of volume and/or mass constraints can be implemented with the introduction of the characterization function to indicate the existence information for each material. Finally, the effectiveness of the proposed methodology is illustrated by three numerical examples.