A method for three-dimensional reconstruction from an irregular ultrasound sampling grid, the method comprising: acquiring ultrasound data representing sample locations in first and second planes associated with scanning from different apertures along a linear array twisted about a longitudinal axis of the linear array such that the sample locations of the first plane are offset from the sample locations of the second plane by a non-integer amount of the distance between adjacent sample locations in the first and/or second planes; determining a grid voxel of a three-dimensional Cartesian grid, the grid voxel between the first and second planes; interpolating a voxel value for the grid voxel; rendering an image from the three-dimensional grid, the image being a function of the voxel value, characterised in that the interpolating comprises: bilinearly interpolating, by a graphics processing unit, a first sample in the first plane from data of the sample locations of the first plane; bilinearly interpolating, by the graphic processing unit, a second sample in the second plane from data of the sample locations of the second plane; linearly interpolating the voxel value for the grid voxel from the first and second samples; wherein the bilinearly interpolating of the first sample comprises interpolating the first sample to a first position on the first plane and closest to the grid voxel, relative distances from the sample locations of the first plane to the first position weighting the interpolating of the first sample; and wherein the bilinearly interpolating of the second sample comprises interpolating the second sample to a second position on the second plane and closest to the grid voxel, relative distances from the sample locations of the second plane to the second position weighting the interpolating of the second sample. acquiring ultrasound data representing sample locations in first and second planes associated with scanning from different apertures along a linear array twisted about a longitudinal axis of the linear array such that the sample locations of the first plane are offset from the sample locations of the second plane by a non-integer amount of the distance between adjacent sample locations in the first and/or second planes; determining a grid voxel of a three-dimensional Cartesian grid, the grid voxel between the first and second planes; interpolating a voxel value for the grid voxel; rendering an image from the three-dimensional grid, the image being a function of the voxel value, characterised in that the interpolating comprises: bilinearly interpolating, by a graphics processing unit, a first sample in the first plane from data of the sample locations of the first plane; bilinearly interpolating, by the graphic processing unit, a second sample in the second plane from data of the sample locations of the second plane; linearly interpolating the voxel value for the grid voxel from the first and second samples; wherein the bilinearly interpolating of the first sample comprises interpolating the first sample to a first position on the first plane and closest to the grid voxel, relative distances from the sample locations of the first plane to the first position weighting the interpolating of the first sample; and wherein the bilinearly interpolating of the second sample comprises interpolating the second sample to a second position on the second plane and closest to the grid voxel, relative distances from the sample locations of the second plane to the second position weighting the interpolating of the second sample. bilinearly interpolating, by a graphics processing unit, a first sample in the first plane from data of the sample locations of the first plane; bilinearly interpolating, by the graphic processing unit, a second sample in the second plane from data of the sample locations of the second plane; linearly interpolating the voxel value for the grid voxel from the first and second samples; wherein the bilinearly interpolating of the first sample comprises interpolating the first sample to a first position on the first plane and closest to the grid voxel, relative distances from the sample locations of the first plane to the first position weighting the interpolating of the first sample; and wherein the bilinearly interpolating of the second sample comprises interpolating the second sample to a second position on the second plane and closest to the grid voxel, relative distances from the sample locations of the second plane to the second position weighting the interpolating of the second sample.