Using the raster package in R you could apply a pixel-wise regression estimate of NDVI ~ time. Here is an example for a linear model, locally-weighted polynomial regression and regression coefficients.

# Create some example data
r <- raster(nrow=100, ncol=100)
  r[] <- runif(ncell(r),-1,1)
    rt <- stack(r)      
      for(i in 2:26) {
        r <- rt[[1]] 
          r[] <- runif(ncell(r),-1,1)     
            rt <- addLayer(rt, r)  
      }

# Create a time vector to act as x
time <- sort(sample(1:365,nlayers(rt))) 
        
# linear (lm) regression estimate(s) of ndvi ~ time
t.lm.predict <- function(x) {if (is.na(x[1])) {NA} else {predict(lm(x ~ time))}} 
f.pred <- calc(rt, t.lm.predict)
  plot(f.pred)

# locally-weighted polynomial regression of ndvi ~ time
t.lowess <- function(x,...) { if (is.na(x[1])) { NA } else { lowess(x,y,...)$y } } 
f.pred <- calc(rt, t.lowess)
  plot(f.pred)
      
# slope and intercept of ndvi ~ time
t.lm.coef <- function(x) {
  if (is.na(x[1])) { NA } else { lm(x ~ time)$coefficients }
  }
f.coef <- calc(rt, t.lm.coef)

Using the raster package in R you could apply a pixel-wise regression estimate of NDVI ~ time. Here is an example for a linear model, locally-weighted polynomial regression and regression coefficients.

library(raster)

# Create some example data
r <- raster(nrow=100, ncol=100)
  r[] <- runif(ncell(r),-1,1)
    rt <- stack(r)      
      for(i in 2:26) {
        r <- rt[[1]] 
          r[] <- runif(ncell(r),-1,1)     
            rt <- addLayer(rt, r)  
      }

# Create a time vector to act as x
time <- sort(sample(1:365,nlayers(rt))) 
        
# linear (lm) regression estimate(s) of ndvi ~ time
t.lm.predict <- function(x) {if (is.na(x[1])) {NA} else {predict(lm(x ~ time))}} 
f.pred <- calc(rt, t.lm.predict)
  plot(f.pred)

# locally-weighted polynomial regression of ndvi ~ time
t.lowess <- function(x,...) { if (is.na(x[1])) { NA } else { lowess(x,y,...)$y } } 
f.pred <- calc(rt, t.lowess)
  plot(f.pred)
      
# slope and intercept of ndvi ~ time
t.lm.coef <- function(x) {
  if (is.na(x[1])) { NA } else { lm(x ~ time)$coefficients }
  }
f.coef <- calc(rt, t.lm.coef)

Using the raster package in R you could apply a pixel-wise regression estimate of NDVI ~ time. Here is an example.

require(raster)

# Create some example data
r <- raster(nrow=100, ncol=100)
  r[] <- runif(ncell(r),-1,1)
    rt <- stack(r)      
      for(i in 2:26) {
        r <- rt[[1]] 
          r[] <- runif(ncell(r),-1,1)     
            rt <- addLayer(rt, r)  
      }

# Create a time vector to act as x
time <- 1:nlayers(rt)   
        
# regression estimate(s) of ndvi ~ time
t.lm.predict <- function(x) {
  if (is.na(x[1])) {
    NA
  } else {
    predict( lm(x ~ time) )
  }
 } 
f.pred <- calc(rt, t.lm.predict)
par(mfrow=c(2,2))
  plot(rt[[1]])
    plot(f.pred[[1]])
       plot(rt[[26]])
         plot(f.pred[[26]])

# slope and intercept of ndvi ~ time
t.lm.coef <- function(x) {
  if (is.na(x[1])) {
    NA
  } else {
    m = lm(x ~ time)
  m$coefficients
  }
}
f.coef <- calc(rt, t.lm.coef)

Using the raster package in R you could apply a pixel-wise regression estimate of NDVI ~ time. Here is an example for a linear model, locally-weighted polynomial regression and regression coefficients.

# Create some example data
r <- raster(nrow=100, ncol=100)
  r[] <- runif(ncell(r),-1,1)
    rt <- stack(r)      
      for(i in 2:26) {
        r <- rt[[1]] 
          r[] <- runif(ncell(r),-1,1)     
            rt <- addLayer(rt, r)  
      }

# Create a time vector to act as x
time <- sort(sample(1:365,nlayers(rt))) 
        
# linear (lm) regression estimate(s) of ndvi ~ time
t.lm.predict <- function(x) {if (is.na(x[1])) {NA} else {predict(lm(x ~ time))}} 
f.pred <- calc(rt, t.lm.predict)
  plot(f.pred) 

# locally-weighted polynomial regression of ndvi ~ time
t.lowess <- function(x,...) { if (is.na(x[1])) { NA } else { lowess(x,y,...)$y } } 
f.pred <- calc(rt, t.lowess)
  plot(f.pred)
      
# slope and intercept of ndvi ~ time
t.lm.coef <- function(x) {
  if (is.na(x[1])) { NA } else { lm(x ~ time)$coefficients }
  }
f.coef <- calc(rt, t.lm.coef)