A measurement of the production cross section of top quark pairs in association with two b jets ($\mathrm{t\bar{t}}\mathrm{b\bar{b}}$) is presented using data collected in proton-proton collisions at $\sqrt{s} =$ 13 TeV by the CMS detector at the LHC corresponding to an integrated luminosity of 35.9 fb$^{-1}$. The cross section is measured in the all-jet decay channel of the top quark pair by selecting events containing at least eight jets, of which at least two are identified as originating from the hadronization of b quarks. A combination of multivariate analysis techniques is used to reduce the large background from multijet events not containing a top quark pair, and to help discriminate between jets originating from top quark decays and other additional jets. The cross section is determined for the total phase space to be 5.5 $\pm$ 0.3 (stat) ${}^{+1.6}_{-1.3}$ (syst) pb and also measured for two fiducial $\mathrm{t\bar{t}}\mathrm{b\bar{b}}$ definitions. The measured cross sections are found to be larger than theoretical predictions by a factor of 1.5-2.4, corresponding to 1-2 standard deviations.