However, at a spatial resolution of about 1 km, the pixels frequently include two or more types of crops and subsequently provide mixed VIs profiles, which might not be the best to fully characterise the spatial patterns of phenological variation relevant to a crop's important growth stages. Due to this, LSP has entered a new era with the availability of data from new missions such as Landsat 8 (L8) and Sentinel-2 (S2). Various studies have generated LSP metrics from the L8 satellite at 30 m resolution [21; 22] but typically lack the temporal and sufficient spatial resolution to capture continuous crop development and reflect the local spatial heterogeneity of crop information. In this way, the Sentinel-2 mission can partially overcome this constraint by offering high spatial resolution of 10 m and a temporal resolution of either 10 days with one satellite or 5 days while using two-satellite constellation [23]. Although S2 boasts various enhanced design features compared to other low- to medium-spatial resolution satellites, in regions with persistent rain and clouds, these optical sensors are naturally impacted by cloud cover, aerosols, and other atmospheric contaminants [24]. This results in undesirable gaps and noise in current VIs products, which pose obstacles for subsequent high-level applications such as LSP analysis and extraction [25]. As a result, one of the essential processes that often precedes any LSP analysis is the construction of high-quality VIs time series using various filtering algorithms such as a Maximum-Value Composite (MVC) [26], the Savitzky-Golay (SG) filter [27; 28], Asymmetric Gaussian (AG) function fitting [25], the best index slope extraction technique (BISE) [29], Double Logistic (DL) function fitting [30], HANTS [31; 32], and the Whittaker smoother [33]. By smoothing the raw times series of VIs and producing high-quality VIs data, LSP detection can be effectively managed. However, accurate extraction of phenology, particularly for crops, is based on identifying valid crop phenological cycles [34]. Despite its theoretical simplicity, this task proves challenging in practice for several reasons. First, while natural vegetation shows a significant potential for periodicity in the variation of their \(NDVI\) seasonal pattern over the years, this is not true when aiming for a crop-specific analysis because multiple crop types with different phenological patterns are frequently planted in rotation on the same plot, resulting in the ability of a selected \(NDVI\) time series to demonstrate the feature of multiple peaks and cycles in a single year [35]. Second, the cultivation time period of a crop type is not always static but dynamically changeable from time to time in response to several variables such as climatic and environmental factors, crop varieties, and interannual variation of planting time, in addition to regional differences. Due to all of these factors, prior knowledge about crop cycle intervals may not be representative of the current crop being grown in a particular pixel; thus, the need arises to develop a robust pixel-based approach that can automatically and accurately determine the crop's growth cycle for each calendar year from the processed time series of VIs data. Afterwards, the LSP parameters can be estimated. Currently, a variety of methods have been developed for the retrieval of phenological parameters from VIs time series. The main idea for most of these methods is to identify the turning points in the time series of satellite-derived VIs. However, they can be categorized into four broad categories: fixed and dynamic threshold, curvature and inflection, trend, and priori curve-based approaches [15]. Among these methods, the threshold-based approach has gained researchers' attention and has become the most widely used approach to estimate both SOS and EOS metrics due to its simplicity and efficiency in generating results within a reasonable time frame compared with other methods [15; 36].