There are a few methods that do what you want, which is to allow functional forms to be flexible. Probably the best one for your case here however is the additive model (or generalized additive model if your response isn't continuous).
The AM has the form
$$
y = \alpha + X'\beta + \displaystyle\sum_m f_m(Z_m) + \epsilon
$$
where $X$ are variables that you represent linearly as in a parametric model, and the $Z$ are terms that you allow to be flexible.
AMs don't remove all functional form assumptions. You still have that your terms have an additive effect on the response, which might not be true (interactions can help with some violations of this). And when they are fit the usual way, they assume that the functions $f_m$ are smooth (so, like, no sharp thresholds).
They are usually fit with a method similar to what you call "segmented linear regression", but modified. One uses way, way more breakpoints than one thinks one needs, and then lets the data tell them how to "penalize" the coefficients on those breakpoints, so that the function doesn't get too wiggly (or too straight).
Specifically, if the response is continuous, the parameters ($\theta$ $\equiv$ $\beta$ along with the splines that make up each of the $f_z$) can be estimated by
$$
\theta = \left(W'W+\lambda D\right)^{-1}W'y
$$
where $W$ is $X$ along with the data gotten from the basis expansion (the segmenting) of $Z$. $D$ is a square matrix of dimension $p$ (the number of variables including the spline terms). In a really simple AM, it is zero on al the off-diagonals, and zero for the $X$'s and $Z$'s, but one for all the spline terms. $\lambda$ is a vector with different values in the entries associated with the sets of splined terms for each $Z$. When the entries of $\lambda$ get bigger, the splined terms shrink, and the function $f(Z)$ reduces to the main effect $Z$ (because in this simple example there is no penalty on $Z$ by the way $D$ is constructed).
Choosing $\lambda$ from data is a little more involved, as is much else about AMs and GAMs. Good books include those by Ruppert, Wand and Carrol ("Semiparametric Regression" -- dated, but easy to read), Wood (introduction to GAMs with R), and the chapter in Hastie & Tibshirani's Intro to Statistical Learning.
While I think there are a few R packages that do GAMs, I've had good luck with gam in mgcv.
Try the following as an example of the kind of output you'd get:
library(mgcv)
m = gam(Sepal.Width~s(Sepal.Length)+s(Petal.Length)+s(Petal.Width)+Species,data=iris)
summary(m)
plot(m,pages=1)
All the above said: GAMs can be kind of a rabbit hole. If you don't have the time and inclination to learn all this stuff, see if your needs aren't served by including a few polynomial terms.
