1. A method for medical imaging of an inclusion comprising: imaging a tissue volume with a probe comprising: an ultrasound transducer that is operative to provide an on-site estimation of inclusion size and location; a first emitter and a first detector; the first emitter having light of a wavelength of about 400 to about 900 nanometers; the first detector detecting light of a wavelength of about 400 to about 900 nanometers; a source circuit connected in operational communication to the emitter; a detector circuit connected in operational communication to the detector; and a central processing unit connected to the source circuit and the detector circuit; scanning the tissue volume with light having a wavelength of about 400 to about 900 nanometers; the tissue volume comprising a first layer and a second layer; segmenting the scanned tissue volume into an inclusion region comprising a plurality of first voxels and a background region comprising a plurality of second voxels; the volume of each second voxel being larger than the volume of each first voxel; identifying a tissue layer thickness and an approximate tilting angle between the tissue layer and the probe; estimating photon density: Ï†(r,Ï‰)=f(Î¼ a1 , Î¼â€² s1 , Î¼ a2 , Î¼â€² s2 ), where Î¼ a1 , Î¼ a2 , Î¼â€² s1 , and Î¼â€² s2 are absorption and reduced scattering coefficients of first and second layers, respectively; wherein absorption coefficients have the subscript â€œaâ€, scattering coefficients have the subscript â€œsâ€, wherein the subscript â€œ1â€ represents the first layer, the subscript â€œ2â€ represents the second layer, and the superscript (â€²) stands for the reduced scattering coefficient; applying a conjugate gradient method to estimate fitted optical properties of the first layer and the second layer; minimizing a least-square objective function given as minâˆ¥Ï†(r,Ï‰)âˆ’Ï† 0 (r,Ï‰)âˆ¥ 2 ; wherein a forward Jacobian weight matrix Wij = [ âˆ‚ Ï• ij âˆ‚ Î¼ aj , âˆ‚ Ï• ij âˆ‚ D j ] , that relates the photon density perturbation at detector i and imaging voxel j with absorption coefficient change Î”Î¼ aj and diffusion coefficient change Î”D j , is calculated by using the bulk optical properties simulated from the two-layer layer model and given in equation (1) as: [ W ij ] = [ âˆ‚ Ï• 11 âˆ‚ Î¼ a â¢ â¢ 1 â€¦ âˆ‚ Ï• 1 â¢ L âˆ‚ Î¼ aL âˆ‚ Ï• 11 âˆ‚ D 1 â€¦ âˆ‚ Ï• 1 â¢ â¢ L âˆ‚ D L âˆ‚ Ï• 21 âˆ‚ Î¼ a â¢ â¢ 1 â€¦ âˆ‚ Ï• 2 â¢ â¢ L âˆ‚ Î¼ aL âˆ‚ Ï• 21 âˆ‚ D 1 â€¦ âˆ‚ Ï• 2 â¢ L âˆ‚ D L â‹® â‹± â‹® â‹® â‹± â‹® âˆ‚ Ï• M â¢ â¢ 1 âˆ‚ Î¼ a â¢ â¢ 1 â€¦ âˆ‚ Ï• ML âˆ‚ Î¼ aL âˆ‚ Ï• M â¢ â¢ 1 âˆ‚ D 1 â€¦ âˆ‚ Ï• ML âˆ‚ D L ] , ( 1 ) where M is the total number of detector readings; L is the total number of imaging voxels and Ï† 0 (r,Ï‰) is the photon density of the first layer.