This paper presents a two-dimensional (2D) fitting algorithm, GALFIT, designed to extract structural components from galaxy images, particularly for well-resolved, nearby galaxies observed with the Hubble Space Telescope (HST). The algorithm improves upon previous techniques by allowing simultaneous fitting of galaxies with an arbitrary number of components and optimizing computation speed for large images. It uses 2D models such as the Nuker law, Sérsic (de Vaucouleurs) profile, exponential disk, and Gaussian or Moffat functions. The algorithm also generalizes azimuthal shapes as generalized ellipses to fit disk-like and boxy components. The algorithm is useful for standard modeling of galaxy profiles, extracting bars, stellar disks, double nuclei, and compact nuclear sources, and measuring absolute dust extinction or surface brightness fluctuations after removing the galaxy model. The study shows that even simple-looking galaxies generally require at least three components for accurate modeling, rather than the one or two components often used. Many galaxies with complex isophotes, ellipticity changes, and position-angle twists can be accurately modeled in 2D. The algorithm is tested on seven case studies, including regular and barred spiral galaxies, highly disk-like lenticular galaxies, and elliptical galaxies with varying complexities. A useful extension is accurately extracting nuclear point sources in galaxies. The algorithm is compared with 1D extraction techniques on simulated images with different nuclear slope cuspiness, and its application to nearby galaxies with weak nuclei is illustrated. The algorithm uses convolution to account for telescope and atmospheric seeing, and it can be turned off if not needed. It minimizes the reduced chi-squared statistic, defined as the sum of squared residuals divided by the degrees of freedom. The algorithm uses a downhill gradient method for parameter optimization and can handle a wide range of galaxy profiles, including Sérsic, exponential disk, Nuker, Gaussian, and Moffat profiles. The algorithm is implemented in C and can be used with FITS image files. It allows for the creation of point sources by approximating a delta function with a Gaussian or Moffat function. The algorithm is also used to create bars by using a Sérsic profile with a flat inner and steep outer profile and a boxy shape. The algorithm is computationally intensive, but it can be optimized by convolving only the area most affected by seeing. The algorithm is tested on a variety of galaxy images, and it is shown to accurately model galaxy profiles with a wide range of parameters. The algorithm is also used to estimate parameter uncertainties by analyzing the covariance matrix and using Monte Carlo simulations to search the parameter space. The algorithm is used to decompose galaxies into components, and the number and types of components are determined iteratively. The algorithm is also used to test for degeneracies in galaxy decompositions, and it is shown that galaxy light profiles can generally be modeled accurately with three to five components. The algorithm is used to study the significance of components in galaxy decomThis paper presents a two-dimensional (2D) fitting algorithm, GALFIT, designed to extract structural components from galaxy images, particularly for well-resolved, nearby galaxies observed with the Hubble Space Telescope (HST). The algorithm improves upon previous techniques by allowing simultaneous fitting of galaxies with an arbitrary number of components and optimizing computation speed for large images. It uses 2D models such as the Nuker law, Sérsic (de Vaucouleurs) profile, exponential disk, and Gaussian or Moffat functions. The algorithm also generalizes azimuthal shapes as generalized ellipses to fit disk-like and boxy components. The algorithm is useful for standard modeling of galaxy profiles, extracting bars, stellar disks, double nuclei, and compact nuclear sources, and measuring absolute dust extinction or surface brightness fluctuations after removing the galaxy model. The study shows that even simple-looking galaxies generally require at least three components for accurate modeling, rather than the one or two components often used. Many galaxies with complex isophotes, ellipticity changes, and position-angle twists can be accurately modeled in 2D. The algorithm is tested on seven case studies, including regular and barred spiral galaxies, highly disk-like lenticular galaxies, and elliptical galaxies with varying complexities. A useful extension is accurately extracting nuclear point sources in galaxies. The algorithm is compared with 1D extraction techniques on simulated images with different nuclear slope cuspiness, and its application to nearby galaxies with weak nuclei is illustrated. The algorithm uses convolution to account for telescope and atmospheric seeing, and it can be turned off if not needed. It minimizes the reduced chi-squared statistic, defined as the sum of squared residuals divided by the degrees of freedom. The algorithm uses a downhill gradient method for parameter optimization and can handle a wide range of galaxy profiles, including Sérsic, exponential disk, Nuker, Gaussian, and Moffat profiles. The algorithm is implemented in C and can be used with FITS image files. It allows for the creation of point sources by approximating a delta function with a Gaussian or Moffat function. The algorithm is also used to create bars by using a Sérsic profile with a flat inner and steep outer profile and a boxy shape. The algorithm is computationally intensive, but it can be optimized by convolving only the area most affected by seeing. The algorithm is tested on a variety of galaxy images, and it is shown to accurately model galaxy profiles with a wide range of parameters. The algorithm is also used to estimate parameter uncertainties by analyzing the covariance matrix and using Monte Carlo simulations to search the parameter space. The algorithm is used to decompose galaxies into components, and the number and types of components are determined iteratively. The algorithm is also used to test for degeneracies in galaxy decompositions, and it is shown that galaxy light profiles can generally be modeled accurately with three to five components. The algorithm is used to study the significance of components in galaxy decom