diff --git a/build_win.bat b/build_win.bat
index dac590811..600d6d6e6 100644
--- a/build_win.bat
+++ b/build_win.bat
@@ -2,67 +2,79 @@
 @IF "%PS_ECHO_ON%" NEQ "" (echo on) ELSE (echo off)
 @GOTO :MAIN
 :HELP
+@ECHO.
 @ECHO Performs initial build or rebuild of the app (build) and deps (build/deps).
 @ECHO Default options are determined from build directories and system state.
 @ECHO.
 @ECHO Usage: build_win [-ARCH ^<arch^>] [-CONFIG ^<config^>] [-DESTDIR ^<directory^>]
 @ECHO                  [-STEPS ^<all^|all-dirty^|app^|app-dirty^|deps^|deps-dirty^>]
+@ECHO                  [-RUN ^<console^|custom^|none^|viewer^|window^>]
 @ECHO.
-@ECHO  -a -ARCH        Target processor architecture
-@ECHO                  Default: %PS_ARCH_HOST%
-@ECHO  -c -CONFIG      MSVC project config
-@ECHO                  Default: %PS_CONFIG_DEFAULT%
-@ECHO  -s -STEPS       Performs only the specified build steps:
-@ECHO                    all - clean and build deps and app
-@ECHO                    all-dirty - build deps and app without cleaning
-@ECHO                    app - build main project/application
-@ECHO                    app-dirty - does not build main project/application
-@ECHO                    deps - clean and build deps
-@ECHO                    deps-dirty - build deps without cleaning
-@ECHO                  Default: %PS_STEPS_DEFAULT%
-@ECHO  -d -DESTDIR     Deps destination directory (ignored on dirty builds)
-@ECHO                  %PS_DESTDIR_DEFAULT_MSG%
+@ECHO  -a -ARCH      Target processor architecture
+@ECHO                Default: %PS_ARCH_HOST%
+@ECHO  -c -CONFIG    MSVC project config
+@ECHO                Default: %PS_CONFIG_DEFAULT%
+@ECHO  -s -STEPS     Performs only the specified build steps:
+@ECHO                  all - clean and build deps and app
+@ECHO                  all-dirty - build deps and app without cleaning
+@ECHO                  app - clean and build main applications
+@ECHO                  app-dirty - build main applications without cleaning
+@ECHO                  deps - clean and build deps
+@ECHO                  deps-dirty - build deps without cleaning
+@ECHO                Default: %PS_STEPS_DEFAULT%
+@ECHO  -r -RUN       Specifies what to perform at the run step:
+@ECHO                  console - run and wait on prusa-slicer-console.exe
+@ECHO                  custom - run and wait on your custom build/%PS_CUSTOM_RUN_FILE%
+@ECHO                  ide - open project in Visual Studio if not open (no wait)
+@ECHO                  none - run step does nothing
+@ECHO                  viewer - run prusa-gcodeviewer.exe (no wait)
+@ECHO                  window - run prusa-slicer.exe (no wait)
+@ECHO                Default: none
+@ECHO  -d -DESTDIR   Deps destination directory
+@ECHO                Warning: Changing destdir path will not delete the old destdir.
+@ECHO                Default: %PS_DESTDIR_DEFAULT_MSG%
 @ECHO.
-@ECHO  Example usage:
-@ECHO                  First build:  build_win -d "c:\src\PrusaSlicer-deps"
-@ECHO                  Deps change:  build_win -s all
-@ECHO                  App rebuild:  build_win
+@ECHO  Examples:
+@ECHO.
+@ECHO  Initial build:           build_win -d "c:\src\PrusaSlicer-deps"
+@ECHO  Build post deps change:  build_win -s all
+@ECHO  App dirty build:         build_win
+@ECHO  App dirty build ^& run:   build_win -r console
+@ECHO  All clean build ^& run:   build_win -s all -r console -d "deps\build\out_deps"
 @ECHO.
 GOTO :END
 
 :MAIN
+REM Script constants
 SET START_TIME=%TIME%
 SET PS_START_DIR=%CD%
+SET PS_SOLUTION_NAME=PrusaSlicer
+SET PS_CHOICE_TIMEOUT=30
+SET PS_CUSTOM_RUN_FILE=custom_run.bat
+SET PS_DEPS_PATH_FILE_NAME=.DEPS_PATH.txt
+SET PS_DEPS_PATH_FILE=%~dp0deps\build\%PS_DEPS_PATH_FILE_NAME%
+SET PS_CONFIG_LIST="Debug;MinSizeRel;Release;RelWithDebInfo"
+
+REM Probe build directories and system state for reasonable default arguments
 pushd %~dp0
-REM Probe build directories and sytem state for reasonable default arguments
 SET PS_CONFIG=RelWithDebInfo
-SET PS_ARCH=%PROCESSOR_ARCHITECTURE%
+SET PS_ARCH=%PROCESSOR_ARCHITECTURE:amd64=x64%
 CALL :TOLOWER PS_ARCH
-SET PS_DEPS_PATH_FILE=%~dp0deps\build\.DEPS_PATH.txt
+SET PS_RUN=none
 SET PS_DESTDIR=
-IF EXIST %PS_DEPS_PATH_FILE% (
-    FOR /F "tokens=* USEBACKQ" %%I IN ("%PS_DEPS_PATH_FILE%") DO SET PS_DESTDIR=%%I
-    IF EXIST build/ALL_BUILD.vcxproj (
-        SET PS_STEPS=app-dirty
-    ) ELSE SET PS_STEPS=app
-) ELSE SET PS_STEPS=all
-SET PS_DESTDIR_CACHED=%PS_DESTDIR%
+CALL :RESOLVE_DESTDIR_CACHE
 
 REM Set up parameters used by help menu
 SET PS_CONFIG_DEFAULT=%PS_CONFIG%
 SET PS_ARCH_HOST=%PS_ARCH%
-SET PS_STEPS_DEFAULT=%PS_STEPS%
-IF "%PS_DESTDIR%" NEQ "" (
-    SET PS_DESTDIR_DEFAULT_MSG=Default: %PS_DESTDIR%
-) ELSE (
-    SET PS_DESTDIR_DEFAULT_MSG=Argument required ^(no default available^)
-)
+(echo " -help /help -h /h -? /? ")| findstr /I /C:" %~1 ">nul && GOTO :HELP
 
 REM Parse arguments
 SET EXIT_STATUS=1
+SET PS_CURRENT_STEP=arguments
 SET PARSER_STATE=
 SET PARSER_FAIL=
-FOR %%I in (%*) DO CALL :PARSE_OPTION "ARCH CONFIG DESTDIR STEPS" PARSER_STATE "%%~I"
+FOR %%I in (%*) DO CALL :PARSE_OPTION "ARCH CONFIG DESTDIR STEPS RUN" PARSER_STATE "%%~I"
 IF "%PARSER_FAIL%" NEQ "" (
     @ECHO ERROR: Invalid switch: %PARSER_FAIL% 1>&2
     GOTO :HELP
@@ -72,46 +84,71 @@ IF "%PARSER_FAIL%" NEQ "" (
 )
 
 REM Validate arguments
-SET PS_STEPS_SAVED=%PS_STEPS%
-CALL :PARSE_OPTION_NAME "all all-dirty deps-dirty deps app-dirty app" PS_STEPS -%PS_STEPS%
-IF "%PS_STEPS%" EQU "" (
-    @ECHO ERROR: Invalid parameter: steps=%PS_STEPS_SAVED% 1>&2
-    GOTO :HELP
-) ELSE SET PS_STEPS_SAVED=
-(echo %PS_STEPS%)| findstr /I /C:"dirty">nul && SET PS_STEPS_DIRTY=1
-CALL :TOLOWER PS_STEPS
+SET PS_ASK_TO_CONTINUE=
 CALL :TOLOWER PS_ARCH
-CALL :CANONICALIZE_PATH PS_DESTDIR "%PS_START_DIR%"
+SET PS_ARCH=%PS_ARCH:amd64=x64%
+CALL :PARSE_OPTION_VALUE %PS_CONFIG_LIST:;= % PS_CONFIG
+IF "%PS_CONFIG%" EQU "" GOTO :HELP
+REM RESOLVE_DESTDIR_CACHE must go after PS_ARCH and PS_CONFIG, but before PS STEPS
+CALL :RESOLVE_DESTDIR_CACHE
+IF "%PS_STEPS%" EQU "" SET PS_STEPS=%PS_STEPS_DEFAULT%
+CALL :PARSE_OPTION_VALUE "all all-dirty deps-dirty deps app-dirty app app-cmake" PS_STEPS
+IF "%PS_STEPS%" EQU "" GOTO :HELP
+(echo %PS_STEPS%)| findstr /I /C:"dirty">nul && SET PS_STEPS_DIRTY=1 || SET PS_STEPS_DIRTY=
+IF "%PS_STEPS%" EQU "app-cmake" SET PS_STEPS_DIRTY=1
+IF "%PS_DESTDIR%" EQU "" SET PS_DESTDIR=%PS_DESTDIR_CACHED%
 IF "%PS_DESTDIR%" EQU "" (
-    IF "%PS_STEPS_DIRTY%" EQU "" (
-        @ECHO ERROR: Parameter required: destdir 1>&2
-        GOTO :HELP
-    )
-) ELSE IF "%PS_DESTDIR%" NEQ "%PS_DESTDIR_CACHED%" (
-    IF "%PS_STEPS_DIRTY%" NEQ "" (
-        @ECHO WARNING: Parameter ignored: destdir
-    ) ELSE (echo "all deps")| findstr /I /C:"%PS_STEPS%">nul || (
-        @ECHO WARNING: Conflict with cached parameter: 1>&2
-        @ECHO WARNING:  -destdir=%PS_DESTDIR% 1>&2
-        @ECHO WARNING:    cached=%PS_DESTDIR_CACHED% 1>&2
+    @ECHO ERROR: Parameter required: -DESTDIR 1>&2
+    GOTO :HELP
+)
+CALL :CANONICALIZE_PATH PS_DESTDIR "%PS_START_DIR%"
+IF "%PS_DESTDIR%" NEQ "%PS_DESTDIR_CACHED%" (
+    (echo "all deps all-dirty deps-dirty")| findstr /I /C:"%PS_STEPS%">nul || (
+        IF EXIST "%PS_DESTDIR%" (
+            @ECHO WARNING: DESTDIR does not match cache: 1>&2
+            @ECHO WARNING:  new: %PS_DESTDIR% 1>&2
+            @ECHO WARNING:  old: %PS_DESTDIR_CACHED% 1>&2
+            SET PS_ASK_TO_CONTINUE=1
+        ) ELSE (
+            @ECHO ERROR: Invalid parameter: DESTDIR=%PS_DESTDIR% 1>&2
+            GOTO :HELP
+        )
     )
 )
 SET PS_DESTDIR_DEFAULT_MSG=
+CALL :PARSE_OPTION_VALUE "console custom ide none viewer window" PS_RUN
+IF "%PS_RUN%" EQU "" GOTO :HELP
+IF "%PS_RUN%" NEQ "none" IF "%PS_STEPS:~0,4%" EQU "deps" (
+    @ECHO ERROR: RUN=%PS_RUN% specified with STEPS=%PS_STEPS%
+    @ECHO ERROR: RUN=none is the only valid option for STEPS "deps" or "deps-dirty"
+    GOTO :HELP
+)
+REM Give the user a chance to cancel if we found something odd.
+IF "%PS_ASK_TO_CONTINUE%" EQU "" GOTO :BUILD_ENV
+@ECHO.
+@ECHO Unexpected parameters detected. Build paused for %PS_CHOICE_TIMEOUT% seconds.
+choice /T %PS_CHOICE_TIMEOUT% /C YN /D N /M "Continue"
+IF %ERRORLEVEL% NEQ 1 GOTO :HELP
 
 REM Set up MSVC environment
+:BUILD_ENV
 SET EXIT_STATUS=2
+SET PS_CURRENT_STEP=environment
 @ECHO **********************************************************************
 @ECHO ** Build Config: %PS_CONFIG%
 @ECHO ** Target Arch:  %PS_ARCH%
 @ECHO ** Build Steps:  %PS_STEPS%
+@ECHO ** Run App:      %PS_RUN%
 @ECHO ** Deps path:    %PS_DESTDIR%
 @ECHO ** Using Microsoft Visual Studio installation found at:
-SET VSWHERE="%ProgramFiles(x86)%\Microsoft Visual Studio\Installer\vswhere.exe"
-IF NOT EXIST %VSWHERE% SET VSWHERE="%ProgramFiles%\Microsoft Visual Studio\Installer\vswhere.exe"
-FOR /F "tokens=* USEBACKQ" %%I IN (`%VSWHERE% -nologo -property installationPath`) DO SET MSVC_DIR=%%I
+SET VSWHERE=%ProgramFiles(x86)%\Microsoft Visual Studio\Installer\vswhere.exe
+IF NOT EXIST "%VSWHERE%" SET VSWHERE=%ProgramFiles%\Microsoft Visual Studio\Installer\vswhere.exe
+FOR /F "tokens=* USEBACKQ" %%I IN (`"%VSWHERE%" -nologo -property installationPath`) DO SET MSVC_DIR=%%I
 @ECHO **  %MSVC_DIR%
 CALL "%MSVC_DIR%\Common7\Tools\vsdevcmd.bat" -arch=%PS_ARCH% -host_arch=%PS_ARCH_HOST% -app_platform=Desktop
-IF "%ERRORLEVEL%" NEQ "0" GOTO :END
+IF %ERRORLEVEL% NEQ 0 GOTO :END
+REM Need to reset the echo state after vsdevcmd.bat clobbers it.
+@IF "%PS_ECHO_ON%" NEQ "" (echo on) ELSE (echo off)
 IF "%PS_DRY_RUN_ONLY%" NEQ "" (
     @ECHO Script terminated early because PS_DRY_RUN_ONLY is set. 1>&2
     GOTO :END
@@ -121,36 +158,132 @@ IF /I "%PS_STEPS:~0,3%" EQU "app" GOTO :BUILD_APP
 REM Build deps
 :BUILD_DEPS
 SET EXIT_STATUS=3
-IF "%PS_STEPS_DIRTY%" EQU "" CALL :MAKE_OR_CLEAN_DIRECTORY deps\build
+SET PS_CURRENT_STEP=deps
+IF "%PS_STEPS_DIRTY%" EQU "" CALL :MAKE_OR_CLEAN_DIRECTORY deps\build "%PS_DEPS_PATH_FILE_NAME%"
 cd deps\build || GOTO :END
-IF "%PS_STEPS_DIRTY%" EQU "" cmake.exe .. -DDESTDIR="%PS_DESTDIR%" || GOTO :END
+cmake.exe .. -DDESTDIR="%PS_DESTDIR%" || GOTO :END
 (echo %PS_DESTDIR%)> "%PS_DEPS_PATH_FILE%"
 msbuild /m ALL_BUILD.vcxproj /p:Configuration=%PS_CONFIG% || GOTO :END
 cd ..\..
-IF /I "%PS_STEPS:~0,4%" EQU "deps" GOTO :PROLOGUE
+IF /I "%PS_STEPS:~0,4%" EQU "deps" GOTO :RUN_APP
 
 REM Build app
 :BUILD_APP
 SET EXIT_STATUS=4
-IF "%PS_STEPS_DIRTY%" EQU "" CALL :MAKE_OR_CLEAN_DIRECTORY build
+SET PS_CURRENT_STEP=app
+IF "%PS_STEPS_DIRTY%" EQU "" CALL :MAKE_OR_CLEAN_DIRECTORY build "%PS_CUSTOM_RUN_FILE%"
 cd build || GOTO :END
-IF "%PS_STEPS_DIRTY%" EQU "" cmake.exe .. -DCMAKE_PREFIX_PATH="%PS_DESTDIR%\usr\local" || GOTO :END
-msbuild /m ALL_BUILD.vcxproj /p:Configuration=%PS_CONFIG% || GOTO :END
+REM Make sure we have a custom batch file skeleton for the run stage
+set PS_CUSTOM_BAT=%PS_CUSTOM_RUN_FILE%
+CALL :CANONICALIZE_PATH PS_CUSTOM_BAT
+IF NOT EXIST %PS_CUSTOM_BAT% CALL :WRITE_CUSTOM_SCRIPT_SKELETON %PS_CUSTOM_BAT%
+SET PS_PROJECT_IS_OPEN=
+FOR /F "tokens=2 delims=," %%I in (
+    'tasklist /V /FI "IMAGENAME eq devenv.exe " /NH /FO CSV ^| find "%PS_SOLUTION_NAME%"'
+) do SET PS_PROJECT_IS_OPEN=%%~I
+cmake.exe .. -DCMAKE_PREFIX_PATH="%PS_DESTDIR%\usr\local" -DCMAKE_CONFIGURATION_TYPES=%PS_CONFIG_LIST% || GOTO :END
+REM Skip the build step if we're using the undocumented app-cmake to regenerate the full config from inside devenv
+IF "%PS_STEPS%" NEQ "app-cmake" msbuild /m ALL_BUILD.vcxproj /p:Configuration=%PS_CONFIG% || GOTO :END
+(echo %PS_DESTDIR%)> "%PS_DEPS_PATH_FILE_FOR_CONFIG%"
+
+REM Run app
+:RUN_APP
+REM All build steps complete.
+CALL :DIFF_TIME ELAPSED_TIME %START_TIME% %TIME%
+IF "%PS_CURRENT_STEP%" NEQ "arguments" (
+    @ECHO.
+    @ECHO Total Build Time Elapsed %ELAPSED_TIME%
+)
+SET EXIT_STATUS=5
+SET PS_CURRENT_STEP=run
+cd src\%PS_CONFIG% || GOTO :END
+IF "%PS_RUN%" EQU "none" GOTO :PROLOGUE
+SET PS_PROJECT_IS_OPEN=
+FOR /F "tokens=2 delims=," %%I in (
+    'tasklist /V /FI "IMAGENAME eq devenv.exe " /NH /FO CSV ^| find "%PS_SOLUTION_NAME%"'
+) do SET PS_PROJECT_IS_OPEN=%%~I
+@ECHO.
+@ECHO Running %PS_RUN% application...
+@REM icacls below is just a hack for file-not-found error handling
+IF "%PS_RUN%" EQU "console" (
+    icacls prusa-slicer-console.exe >nul || GOTO :END
+    start /wait /b prusa-slicer-console.exe
+) ELSE IF "%PS_RUN%" EQU "window" (
+    icacls prusa-slicer.exe >nul || GOTO :END
+    start prusa-slicer.exe
+) ELSE IF "%PS_RUN%" EQU "viewer" (
+    icacls prusa-gcodeviewer.exe >nul || GOTO :END
+    start prusa-gcodeviewer.exe
+) ELSE IF "%PS_RUN%" EQU "custom" (
+    icacls %PS_CUSTOM_BAT% >nul || GOTO :END
+    CALL %PS_CUSTOM_BAT%
+) ELSE IF "%PS_RUN%" EQU "ide" (
+    IF "%PS_PROJECT_IS_OPEN%" NEQ "" (
+        @ECHO WARNING: Solution is already open in Visual Studio. Skipping ide run step. 1>&2
+    ) ELSE (
+        @ECHO Preparing to run Visual Studio...
+        cd ..\.. || GOTO :END
+        REM This hack generates a single config for MSVS, guaranteeing it gets set as the active config.
+        cmake.exe .. -DCMAKE_PREFIX_PATH="%PS_DESTDIR%\usr\local" -DCMAKE_CONFIGURATION_TYPES=%PS_CONFIG% > nul 2> nul || GOTO :END
+        REM Now launch devenv with the single config (setting it active) and a /command switch to re-run cmake and generate the full config list
+        start devenv.exe %PS_SOLUTION_NAME%.sln /command ^"shell /o ^^^"%~f0^^^" -d ^^^"%PS_DESTDIR%^^^" -c %PS_CONFIG% -a %PS_ARCH% -r none -s app-cmake^"
+        REM If devenv fails to launch just directly regenerate the full config list.
+        IF %ERRORLEVEL% NEQ 0 (
+            cmake.exe .. -DCMAKE_PREFIX_PATH="%PS_DESTDIR%\usr\local" -DCMAKE_CONFIGURATION_TYPES=%PS_CONFIG_LIST% 2> nul 1> nul || GOTO :END
+        )
+    )
+)
+
+@REM **********    DON'T ADD ANY CODE BETWEEN THESE TWO SECTIONS    **********
+@REM RUN_APP may hand off control, so let exit codes fall through to PROLOGUE.
 
 :PROLOGUE
 SET EXIT_STATUS=%ERRORLEVEL%
 :END
 @IF "%PS_ECHO_ON%%PS_DRY_RUN_ONLY%" NEQ "" (
-    @ECHO Script Parameters:
+    @ECHO **********************************************************************
+    @ECHO ** Script Parameters:
+    @ECHO **********************************************************************
     @SET PS_
 )
-@ECHO Script started at %START_TIME% and completed at %TIME%.
+IF "%EXIT_STATUS%" NEQ "0" (
+    IF "%PS_CURRENT_STEP%" NEQ "arguments" (
+        @ECHO.
+        @ECHO ERROR: *** Build process failed at %PS_CURRENT_STEP% step. *** 1>&2
+    )
+) ELSE (
+    @ECHO All steps completed successfully.
+)
 popd
 exit /B %EXIT_STATUS%
 
 GOTO :EOF
 REM Functions and stubs start here.
 
+:RESOLVE_DESTDIR_CACHE
+@REM Resolves all DESTDIR cache values and sets PS_STEPS_DEFAULT
+@REM Note: This just sets global variableq, so it doesn't use setlocal.
+SET PS_DEPS_PATH_FILE_FOR_CONFIG=%~dp0build\%PS_ARCH%\%PS_CONFIG%\%PS_DEPS_PATH_FILE_NAME%
+CALL :CANONICALIZE_PATH PS_DEPS_PATH_FILE_FOR_CONFIG
+IF EXIST "%PS_DEPS_PATH_FILE_FOR_CONFIG%" (
+    FOR /F "tokens=* USEBACKQ" %%I IN ("%PS_DEPS_PATH_FILE_FOR_CONFIG%") DO (
+        SET PS_DESTDIR_CACHED=%%I
+        SET PS_DESTDIR_DEFAULT_MSG=%%I
+    )
+    SET PS_STEPS_DEFAULT=app-dirty
+) ELSE IF EXIST "%PS_DEPS_PATH_FILE%" (
+    FOR /F "tokens=* USEBACKQ" %%I IN ("%PS_DEPS_PATH_FILE%") DO (
+        SET PS_DESTDIR_CACHED=%%I
+        SET PS_DESTDIR_DEFAULT_MSG=%%I
+    )
+    SET PS_STEPS_DEFAULT=app
+) ELSE (
+    SET PS_DESTDIR_CACHED=
+    SET PS_DESTDIR_DEFAULT_MSG=Cache missing. Argument required.
+    SET PS_STEPS_DEFAULT=all
+)
+GOTO :EOF
+
 :PARSE_OPTION
 @REM Argument parser called for each argument
 @REM %1 - Valid option list
@@ -179,12 +312,27 @@ IF "%LAST_ARG%%ARG_TYPE%" EQU "NAME" SET PARSER_FAIL=%~3
 )
 GOTO :EOF
 
+:PARSE_OPTION_VALUE
+setlocal disableDelayedExpansion
+@REM Parses value and verifies it is within the supplied list
+@REM %1 - Valid option list
+@REM %2 - In/out variable name; unset on error
+CALL SET NAME=%~2
+CALL SET SAVED_VALUE=%%%NAME%%%
+CALL :PARSE_OPTION_NAME %1 %NAME% -%SAVED_VALUE%
+CALL SET NEW_VALUE=%%%NAME%%%
+IF "%NEW_VALUE%" EQU "" (
+    @ECHO ERROR: Invalid parameter: %NAME:~3%=%SAVED_VALUE% 1>&2
+)
+endlocal & SET %NAME%=%NEW_VALUE%
+GOTO :EOF
+
 :PARSE_OPTION_NAME
 @REM Parses an option name
 @REM %1 - Valid option list
 @REM %2 - Out variable name; unset on error
 @REM %3 - Current argument value
-@REM $4 - Boolean indicating single character switches are valid
+@REM %4 - Boolean indicating single character switches are valid
 @REM Note: Delayed expansion safe because ! character is invalid in option name
 setlocal enableDelayedExpansion
 IF "%4" NEQ "" FOR %%I IN (%~1) DO @(
@@ -201,6 +349,7 @@ IF "%4" NEQ "" (
 )
 @(echo %OPTION_NAME%)| findstr /R /C:".[ ][ ]*.">nul && GOTO :PARSE_OPTION_NAME_FAIL
 @(echo  %~1 )| findstr /I /C:" %OPTION_NAME% ">nul || GOTO :PARSE_OPTION_NAME_FAIL
+FOR %%I IN (%~1) DO SET OPTION_NAME=!OPTION_NAME:%%~I=%%~I!
 endlocal & SET %~2=%OPTION_NAME%
 GOTO :EOF
 :PARSE_OPTION_NAME_FAIL
@@ -210,14 +359,28 @@ GOTO :EOF
 :MAKE_OR_CLEAN_DIRECTORY
 @REM Create directory if it doesn't exist or clean it if it does
 @REM %1 - Directory path to clean or create
+@REM %* - Optional list of files/dirs to keep (in the base directory only)
 setlocal disableDelayedExpansion
 IF NOT EXIST "%~1" (
-    ECHO Creating %~1
-    mkdir "%~1" && GOTO :EOF
+    @ECHO Creating %~1
+    mkdir "%~1" && (
+        endlocal
+        GOTO :EOF
+    )
+)
+@ECHO Cleaning %~1 ...
+SET KEEP_LIST=
+:MAKE_OR_CLEAN_DIRECTORY_ARG_LOOP
+IF "%~2" NEQ "" (
+    SET KEEP_LIST=%KEEP_LIST% "%~2"
+    SHIFT /2
+    GOTO :MAKE_OR_CLEAN_DIRECTORY_ARG_LOOP
 )
-ECHO Cleaning %~1 ...
 for /F "usebackq delims=" %%I in (`dir /a /b "%~1"`) do (
-    (rmdir /s /q  "%~1\%%I" 2>nul ) || del /q /f "%~1\%%I")
+    (echo %KEEP_LIST%)| findstr /I /L /C:"\"%%I\"">nul || (
+        rmdir /s /q  "%~1\%%I" 2>nul ) || del /q /f "%~1\%%I"
+)
+endlocal
 GOTO :EOF
 
 :TOLOWER
@@ -239,7 +402,51 @@ CALL :CANONICALIZE_PATH_INNER %1 %%%~1%% %2
 endlocal & SET %~1=%OUTPUT%
 GOTO :EOF
 :CANONICALIZE_PATH_INNER
-if "%~3" NEQ "" (pushd %~3 || GOTO :EOF)
+if "%~3" NEQ "" (pushd %3 || GOTO :EOF)
 SET OUTPUT=%~f2
-if "%~3" NEQ "" popd %~3
+if "%~3" NEQ "" popd
+GOTO :EOF
+
+:DIFF_TIME
+@REM Calculates elapsed time between two timestamps (TIME environment variable format)
+@REM %1 - Output variable
+@REM %2 - Start time
+@REM %3 - End time
+setlocal EnableDelayedExpansion
+set START_ARG=%2
+set END_ARG=%3
+set END=!END_ARG:%TIME:~8,1%=%%100)*100+1!
+set START=!START_ARG:%TIME:~8,1%=%%100)*100+1!
+set /A DIFF=((((10!END:%TIME:~2,1%=%%100)*60+1!%%100)-((((10!START:%TIME:~2,1%=%%100)*60+1!%%100), DIFF-=(DIFF^>^>31)*24*60*60*100
+set /A CC=DIFF%%100+100,DIFF/=100,SS=DIFF%%60+100,DIFF/=60,MM=DIFF%%60+100,HH=DIFF/60+100
+@endlocal & set %1=%HH:~1%%TIME:~2,1%%MM:~1%%TIME:~2,1%%SS:~1%%TIME:~8,1%%CC:~1%
+@GOTO :EOF
+
+:WRITE_CUSTOM_SCRIPT_SKELETON
+@REM Writes the following text to the supplied file
+@REM %1 - Output filename
+setlocal
+@(
+ECHO @ECHO.
+ECHO @ECHO ********************************************************************************
+ECHO @ECHO ** This is a custom run script skeleton.
+ECHO @ECHO ********************************************************************************
+ECHO @ECHO.
+ECHO @ECHO ********************************************************************************
+ECHO @ECHO ** The working directory is:
+ECHO @ECHO ********************************************************************************
+ECHO dir
+ECHO @ECHO.
+ECHO @ECHO ********************************************************************************
+ECHO @ECHO ** The environment is:
+ECHO @ECHO ********************************************************************************
+ECHO set
+ECHO @ECHO.
+ECHO @ECHO ********************************************************************************
+ECHO @ECHO ** Edit or replace this script to run custom steps after a successful build:
+ECHO @ECHO ** %~1
+ECHO @ECHO ********************************************************************************
+ECHO @ECHO.
+) > "%~1"
+endlocal
 GOTO :EOF
diff --git a/src/libslic3r/EdgeGrid.hpp b/src/libslic3r/EdgeGrid.hpp
index 651a6d763..3c9929149 100644
--- a/src/libslic3r/EdgeGrid.hpp
+++ b/src/libslic3r/EdgeGrid.hpp
@@ -20,10 +20,13 @@ public:
 	Contour(const Slic3r::Point *data, size_t size, bool open) : Contour(data, data + size, open) {}
 	Contour(const std::vector<Slic3r::Point> &pts, bool open) : Contour(pts.data(), pts.size(), open) {}
 
-	const Slic3r::Point *begin()  const { return m_begin; }
-	const Slic3r::Point *end()    const { return m_end; }
-	bool                 open()   const { return m_open; }
-	bool                 closed() const { return ! m_open; }
+    const Slic3r::Point *begin()  const { return m_begin; }
+    const Slic3r::Point *end()    const { return m_end; }
+    bool                 open()   const { return m_open; }
+    bool                 closed() const { return !m_open; }
+
+    const Slic3r::Point &front()  const { return *m_begin; }
+    const Slic3r::Point &back()   const { return *(m_end - 1); }
 
 	// Start point of a segment idx.
 	const Slic3r::Point& segment_start(size_t idx) const {
@@ -61,6 +64,23 @@ public:
 
 	size_t               num_segments() const { return this->size() - (m_open ? 1 : 0); }
 
+    Line                 get_segment(size_t idx) const
+    {
+        assert(idx < this->num_segments());
+        return Line(this->segment_start(idx), this->segment_end(idx));
+    }
+
+    Lines                get_segments() const
+    {
+        Lines lines;
+        lines.reserve(this->num_segments());
+        if (this->num_segments() > 2) {
+            for (auto it = this->begin(); it != this->end() - 1; ++it) lines.push_back(Line(*it, *(it + 1)));
+            if (!m_open) lines.push_back(Line(this->back(), this->front()));
+        }
+        return lines;
+    }
+
 private:
 	size_t  			 size() const { return m_end - m_begin; }
 
diff --git a/src/libslic3r/Layer.cpp b/src/libslic3r/Layer.cpp
index 3f2327686..39228516c 100644
--- a/src/libslic3r/Layer.cpp
+++ b/src/libslic3r/Layer.cpp
@@ -4,6 +4,7 @@
 #include "Fill/Fill.hpp"
 #include "ShortestPath.hpp"
 #include "SVG.hpp"
+#include "BoundingBox.hpp"
 
 #include <boost/log/trivial.hpp>
 
@@ -258,4 +259,26 @@ void Layer::export_region_fill_surfaces_to_svg_debug(const char *name) const
     this->export_region_fill_surfaces_to_svg(debug_out_path("Layer-fill_surfaces-%s-%d.svg", name, idx ++).c_str());
 }
 
+BoundingBox get_extents(const LayerRegion &layer_region)
+{
+    BoundingBox bbox;
+    if (!layer_region.slices.surfaces.empty()) {
+        bbox = get_extents(layer_region.slices.surfaces.front());
+        for (auto it = layer_region.slices.surfaces.cbegin() + 1; it != layer_region.slices.surfaces.cend(); ++it)
+            bbox.merge(get_extents(*it));
+    }
+    return bbox;
+}
+
+BoundingBox get_extents(const LayerRegionPtrs &layer_regions)
+{
+    BoundingBox bbox;
+    if (!layer_regions.empty()) {
+        bbox = get_extents(*layer_regions.front());
+        for (auto it = layer_regions.begin() + 1; it != layer_regions.end(); ++it)
+            bbox.merge(get_extents(**it));
+    }
+    return bbox;
+}
+
 }
diff --git a/src/libslic3r/Layer.hpp b/src/libslic3r/Layer.hpp
index 6cabdeb40..e6bf4b4fc 100644
--- a/src/libslic3r/Layer.hpp
+++ b/src/libslic3r/Layer.hpp
@@ -211,6 +211,9 @@ inline std::vector<float> zs_from_layers(const LayerContainer &layers)
     return zs;
 }
 
+extern BoundingBox get_extents(const LayerRegion &layer_region);
+extern BoundingBox get_extents(const LayerRegionPtrs &layer_regions);
+
 }
 
 #endif
diff --git a/src/libslic3r/Model.cpp b/src/libslic3r/Model.cpp
index 5c291a0c6..f6980e5b2 100644
--- a/src/libslic3r/Model.cpp
+++ b/src/libslic3r/Model.cpp
@@ -1913,14 +1913,16 @@ arrangement::ArrangePolygon ModelInstance::get_arrange_polygon() const
 indexed_triangle_set FacetsAnnotation::get_facets(const ModelVolume& mv, EnforcerBlockerType type) const
 {
     TriangleSelector selector(mv.mesh());
-    selector.deserialize(m_data);
+    // Reset of TriangleSelector is done inside TriangleSelector's constructor, so we don't need it to perform it again in deserialize().
+    selector.deserialize(m_data, false);
     return selector.get_facets(type);
 }
 
 indexed_triangle_set FacetsAnnotation::get_facets_strict(const ModelVolume& mv, EnforcerBlockerType type) const
 {
     TriangleSelector selector(mv.mesh());
-    selector.deserialize(m_data);
+    // Reset of TriangleSelector is done inside TriangleSelector's constructor, so we don't need it to perform it again in deserialize().
+    selector.deserialize(m_data, false);
     return selector.get_facets_strict(type);
 }
 
diff --git a/src/libslic3r/MultiMaterialSegmentation.cpp b/src/libslic3r/MultiMaterialSegmentation.cpp
index 9dde9e9a9..65284ffac 100644
--- a/src/libslic3r/MultiMaterialSegmentation.cpp
+++ b/src/libslic3r/MultiMaterialSegmentation.cpp
@@ -12,6 +12,8 @@
 
 #include <boost/log/trivial.hpp>
 #include <tbb/parallel_for.h>
+#include <mutex>
+#include <boost/thread/lock_guard.hpp>
 
 namespace Slic3r {
 struct ColoredLine {
@@ -38,6 +40,10 @@ struct segment_traits<Slic3r::ColoredLine> {
 };
 }
 
+//#define MMU_SEGMENTATION_DEBUG_GRAPH
+//#define MMU_SEGMENTATION_DEBUG_REGIONS
+//#define MMU_SEGMENTATION_DEBUG_INPUT
+
 namespace Slic3r {
 
 // Assumes that is at most same projected_l length or below than projection_l
@@ -74,7 +80,7 @@ struct PaintedLine
 
 struct PaintedLineVisitor
 {
-    PaintedLineVisitor(const EdgeGrid::Grid &grid, std::vector<PaintedLine> &painted_lines, size_t reserve) : grid(grid), painted_lines(painted_lines)
+    PaintedLineVisitor(const EdgeGrid::Grid &grid, std::vector<PaintedLine> &painted_lines, std::mutex &painted_lines_mutex, size_t reserve) : grid(grid), painted_lines(painted_lines), painted_lines_mutex(painted_lines_mutex)
     {
         painted_lines_set.reserve(reserve);
     }
@@ -115,8 +121,11 @@ struct PaintedLineVisitor
                         if ((line_to_test_projected.a - grid_line.a).cast<double>().squaredNorm() > (line_to_test_projected.b - grid_line.a).cast<double>().squaredNorm())
                             line_to_test_projected.reverse();
 
-                        painted_lines.push_back({it_contour_and_segment->first, it_contour_and_segment->second, line_to_test_projected, this->color});
                         painted_lines_set.insert(*it_contour_and_segment);
+                        {
+                            boost::lock_guard<std::mutex> lock(painted_lines_mutex);
+                            painted_lines.push_back({it_contour_and_segment->first, it_contour_and_segment->second, line_to_test_projected, this->color});
+                        }
                     }
                 }
             }
@@ -127,6 +136,7 @@ struct PaintedLineVisitor
 
     const EdgeGrid::Grid                                                                 &grid;
     std::vector<PaintedLine>                                                             &painted_lines;
+    std::mutex                                                                           &painted_lines_mutex;
     Line                                                                                  line_to_test;
     std::unordered_set<std::pair<size_t, size_t>, boost::hash<std::pair<size_t, size_t>>> painted_lines_set;
     int                                                                                   color             = -1;
@@ -136,14 +146,14 @@ struct PaintedLineVisitor
     static inline const double                                                            append_threshold2 = Slic3r::sqr(append_threshold);
 };
 
-static std::vector<ColoredLine> to_colored_lines(const Polygon &polygon, int color)
+static std::vector<ColoredLine> to_colored_lines(const EdgeGrid::Contour &contour, int color)
 {
     std::vector<ColoredLine> lines;
-    if (polygon.points.size() > 2) {
-        lines.reserve(polygon.points.size());
-        for (auto it = polygon.points.begin(); it != polygon.points.end() - 1; ++it)
+    if (contour.num_segments() > 2) {
+        lines.reserve(contour.num_segments());
+        for (auto it = contour.begin(); it != contour.end() - 1; ++it)
             lines.push_back({Line(*it, *(it + 1)), color});
-        lines.push_back({Line(polygon.points.back(), polygon.points.front()), color});
+        lines.push_back({Line(contour.back(), contour.front()), color});
     }
     return lines;
 }
@@ -238,7 +248,9 @@ static std::vector<ColoredLine> colorize_line(const Line &              line_to_
                                               std::vector<PaintedLine> &painted_lines)
 {
     std::vector<PaintedLine> internal_painted;
-    for (size_t line_idx = start_idx; line_idx <= end_idx; ++line_idx) { internal_painted.emplace_back(painted_lines[line_idx]); }
+    for (size_t line_idx = start_idx; line_idx <= end_idx; ++line_idx)
+        internal_painted.emplace_back(painted_lines[line_idx]);
+
     const int                filter_eps_value = scale_(0.1f);
     std::vector<PaintedLine> filtered_lines;
     filtered_lines.emplace_back(internal_painted.front());
@@ -324,18 +336,18 @@ static std::vector<ColoredLine> colorize_line(const Line &              line_to_
             if (line_1.line.length() <= scale_(0.2)) line_1.color = line_0.color;
     }
 
-    std::vector<ColoredLine> colored_lines_simpl;
-    colored_lines_simpl.emplace_back(final_lines.front());
+    std::vector<ColoredLine> colored_lines_simple;
+    colored_lines_simple.emplace_back(final_lines.front());
     for (size_t line_idx = 1; line_idx < final_lines.size(); ++line_idx) {
         const ColoredLine &line_0 = final_lines[line_idx];
 
-        if (colored_lines_simpl.back().color == line_0.color)
-            colored_lines_simpl.back().line.b = line_0.line.b;
+        if (colored_lines_simple.back().color == line_0.color)
+            colored_lines_simple.back().line.b = line_0.line.b;
         else
-            colored_lines_simpl.emplace_back(line_0);
+            colored_lines_simple.emplace_back(line_0);
     }
 
-    final_lines = colored_lines_simpl;
+    final_lines = colored_lines_simple;
 
     if (final_lines.size() > 1) {
         if (final_lines.front().color != final_lines[1].color && final_lines.front().line.length() <= scale_(0.2)) {
@@ -354,13 +366,12 @@ static std::vector<ColoredLine> colorize_line(const Line &              line_to_
     return final_lines;
 }
 
-static std::vector<ColoredLine> colorize_polygon(const Polygon &poly, const size_t start_idx, const size_t end_idx, std::vector<PaintedLine> &painted_lines)
+static std::vector<ColoredLine> colorize_polygon(const EdgeGrid::Contour &contour, const size_t start_idx, const size_t end_idx, std::vector<PaintedLine> &painted_lines)
 {
     std::vector<ColoredLine> new_lines;
-    Lines                    lines = poly.lines();
-
+    new_lines.reserve(end_idx - start_idx + 1);
     for (size_t idx = 0; idx < painted_lines[start_idx].line_idx; ++idx)
-        new_lines.emplace_back(ColoredLine{lines[idx], 0});
+        new_lines.emplace_back(ColoredLine{contour.get_segment(idx), 0});
 
     for (size_t first_idx = start_idx; first_idx <= end_idx; ++first_idx) {
         size_t second_idx = first_idx;
@@ -368,18 +379,18 @@ static std::vector<ColoredLine> colorize_polygon(const Polygon &poly, const size
         --second_idx;
 
         assert(painted_lines[first_idx].line_idx == painted_lines[second_idx].line_idx);
-        std::vector<ColoredLine> lines_c_line = colorize_line(lines[painted_lines[first_idx].line_idx], first_idx, second_idx, painted_lines);
+        std::vector<ColoredLine> lines_c_line = colorize_line(contour.get_segment(painted_lines[first_idx].line_idx), first_idx, second_idx, painted_lines);
         new_lines.insert(new_lines.end(), lines_c_line.begin(), lines_c_line.end());
 
         if (second_idx + 1 <= end_idx)
             for (size_t idx = painted_lines[second_idx].line_idx + 1; idx < painted_lines[second_idx + 1].line_idx; ++idx)
-                new_lines.emplace_back(ColoredLine{lines[idx], 0});
+                new_lines.emplace_back(ColoredLine{contour.get_segment(idx), 0});
 
         first_idx = second_idx;
     }
 
-    for (size_t idx = painted_lines[end_idx].line_idx + 1; idx < poly.size(); ++idx)
-        new_lines.emplace_back(ColoredLine{lines[idx], 0});
+    for (size_t idx = painted_lines[end_idx].line_idx + 1; idx < contour.num_segments(); ++idx)
+        new_lines.emplace_back(ColoredLine{contour.get_segment(idx), 0});
 
     for (size_t line_idx = 2; line_idx < new_lines.size(); ++line_idx) {
         const ColoredLine &line_0 = new_lines[line_idx - 2];
@@ -456,15 +467,16 @@ static std::vector<ColoredLine> colorize_polygon(const Polygon &poly, const size
     return new_lines;
 }
 
-static std::vector<std::vector<ColoredLine>> colorize_polygons(const Polygons &polygons, std::vector<PaintedLine> &painted_lines)
+static std::vector<std::vector<ColoredLine>> colorize_polygons(const std::vector<EdgeGrid::Contour> &contours, std::vector<PaintedLine> &painted_lines)
 {
     const size_t start_idx = 0;
     const size_t end_idx   = painted_lines.size() - 1;
 
     std::vector<std::vector<ColoredLine>> new_polygons;
+    new_polygons.reserve(contours.size());
 
     for (size_t idx = 0; idx < painted_lines[start_idx].contour_idx; ++idx)
-        new_polygons.emplace_back(to_colored_lines(polygons[idx], 0));
+        new_polygons.emplace_back(to_colored_lines(contours[idx], 0));
 
     for (size_t first_idx = start_idx; first_idx <= end_idx; ++first_idx) {
         size_t second_idx = first_idx;
@@ -473,18 +485,17 @@ static std::vector<std::vector<ColoredLine>> colorize_polygons(const Polygons &p
         --second_idx;
 
         assert(painted_lines[first_idx].contour_idx == painted_lines[second_idx].contour_idx);
-        std::vector<ColoredLine> polygon_c = colorize_polygon(polygons[painted_lines[first_idx].contour_idx], first_idx, second_idx, painted_lines);
-        new_polygons.emplace_back(polygon_c);
+        new_polygons.emplace_back(colorize_polygon(contours[painted_lines[first_idx].contour_idx], first_idx, second_idx, painted_lines));
 
         if (second_idx + 1 <= end_idx)
             for (size_t idx = painted_lines[second_idx].contour_idx + 1; idx < painted_lines[second_idx + 1].contour_idx; ++idx)
-                new_polygons.emplace_back(to_colored_lines(polygons[idx], 0));
+                new_polygons.emplace_back(to_colored_lines(contours[idx], 0));
 
         first_idx = second_idx;
     }
 
-    for (size_t idx = painted_lines[end_idx].contour_idx + 1; idx < polygons.size(); ++idx)
-        new_polygons.emplace_back(to_colored_lines(polygons[idx], 0));
+    for (size_t idx = painted_lines[end_idx].contour_idx + 1; idx < contours.size(); ++idx)
+        new_polygons.emplace_back(to_colored_lines(contours[idx], 0));
 
     return new_polygons;
 }
@@ -507,7 +518,6 @@ struct MMU_Graph
         size_t   to_idx;
         int      color;
         ARC_TYPE type;
-        bool     used{false};
 
         bool operator==(const Arc &rhs) const { return (from_idx == rhs.from_idx) && (to_idx == rhs.to_idx) && (color == rhs.color) && (type == rhs.type); }
         bool operator!=(const Arc &rhs) const { return !operator==(rhs); }
@@ -515,15 +525,16 @@ struct MMU_Graph
 
     struct Node
     {
-        Point                     point;
-        std::list<MMU_Graph::Arc> neighbours;
+        Point             point;
+        std::list<size_t> arc_idxs;
 
-        void remove_edge(const size_t to_idx)
+        void remove_edge(const size_t to_idx, MMU_Graph &graph)
         {
-            for (auto arc_it = this->neighbours.begin(); arc_it != this->neighbours.end(); ++arc_it) {
-                if (arc_it->to_idx == to_idx) {
-                    assert(arc_it->type != ARC_TYPE::BORDER);
-                    this->neighbours.erase(arc_it);
+            for (auto arc_it = this->arc_idxs.begin(); arc_it != this->arc_idxs.end(); ++arc_it) {
+                MMU_Graph::Arc &arc = graph.arcs[*arc_it];
+                if (arc.to_idx == to_idx) {
+                    assert(arc.type != ARC_TYPE::BORDER);
+                    this->arc_idxs.erase(arc_it);
                     break;
                 }
             }
@@ -539,8 +550,8 @@ struct MMU_Graph
 
     void remove_edge(const size_t from_idx, const size_t to_idx)
     {
-        nodes[from_idx].remove_edge(to_idx);
-        nodes[to_idx].remove_edge(from_idx);
+        nodes[from_idx].remove_edge(to_idx, *this);
+        nodes[to_idx].remove_edge(from_idx, *this);
     }
 
     [[nodiscard]] size_t get_global_index(const size_t poly_idx, const size_t point_idx) const { return polygon_idx_offset[poly_idx] + point_idx; }
@@ -548,42 +559,55 @@ struct MMU_Graph
     void append_edge(const size_t &from_idx, const size_t &to_idx, int color = -1, ARC_TYPE type = ARC_TYPE::NON_BORDER)
     {
         // Don't append duplicate edges between the same nodes.
-        for (const MMU_Graph::Arc &arc : this->nodes[from_idx].neighbours)
-            if (arc.to_idx == to_idx)
+        for (const size_t &arc_idx : this->nodes[from_idx].arc_idxs)
+            if (arcs[arc_idx].to_idx == to_idx)
                 return;
-        for (const MMU_Graph::Arc &arc : this->nodes[to_idx].neighbours)
-            if (arc.to_idx == to_idx)
+        for (const size_t &arc_idx : this->nodes[to_idx].arc_idxs)
+            if (arcs[arc_idx].to_idx == to_idx)
                 return;
 
-        this->nodes[from_idx].neighbours.push_back({from_idx, to_idx, color, type});
-        this->nodes[to_idx].neighbours.push_back({to_idx, from_idx, color, type});
+        this->nodes[from_idx].arc_idxs.push_back(this->arcs.size());
         this->arcs.push_back({from_idx, to_idx, color, type});
-        this->arcs.push_back({to_idx, from_idx, color, type});
+
+        // Always insert only one directed arc for the input polygons.
+        // Two directed arcs in both directions are inserted if arcs aren't between points of the input polygons.
+        if (type == ARC_TYPE::NON_BORDER) {
+            this->nodes[to_idx].arc_idxs.push_back(this->arcs.size());
+            this->arcs.push_back({to_idx, from_idx, color, type});
+        }
     }
 
-    // Ignoring arcs in the opposite direction
-    MMU_Graph::Arc get_arc(size_t idx) { return this->arcs[idx * 2]; }
+    // It assumes that between points of the input polygons is always only one directed arc,
+    // with the same direction as lines of the input polygon.
+    [[nodiscard]] MMU_Graph::Arc get_border_arc(size_t idx) const {
+        assert(idx < this->all_border_points);
+        return this->arcs[idx];
+    }
 
     [[nodiscard]] size_t nodes_count() const { return this->nodes.size(); }
 
     void remove_nodes_with_one_arc()
     {
         std::queue<size_t> update_queue;
-        for (const MMU_Graph::Node &node : this->nodes)
-            if (node.neighbours.size() == 1) update_queue.emplace(&node - &this->nodes.front());
+        for (const MMU_Graph::Node &node : this->nodes) {
+            size_t node_idx = &node - &this->nodes.front();
+            // Skip nodes that represent points of input polygons.
+            if (node.arc_idxs.size() == 1 && node_idx >= this->all_border_points)
+                update_queue.emplace(&node - &this->nodes.front());
+        }
 
         while (!update_queue.empty()) {
             size_t           node_from_idx = update_queue.front();
             MMU_Graph::Node &node_from     = this->nodes[update_queue.front()];
             update_queue.pop();
-            if (node_from.neighbours.empty())
+            if (node_from.arc_idxs.empty())
                 continue;
 
-            assert(node_from.neighbours.size() == 1);
-            size_t           node_to_idx = node_from.neighbours.front().to_idx;
+            assert(node_from.arc_idxs.size() == 1);
+            size_t           node_to_idx = arcs[node_from.arc_idxs.front()].to_idx;
             MMU_Graph::Node &node_to     = this->nodes[node_to_idx];
             this->remove_edge(node_from_idx, node_to_idx);
-            if (node_to.neighbours.size() == 1)
+            if (node_to.arc_idxs.size() == 1 && node_to_idx >= this->all_border_points)
                 update_queue.emplace(node_to_idx);
         }
     }
@@ -660,17 +684,17 @@ struct MMU_Graph
             vertex.color(-1);
             Point vertex_point = mk_point(vertex);
 
-            const Point &first_point  = this->nodes[this->get_arc(vertex.incident_edge()->cell()->source_index()).from_idx].point;
-            const Point &second_point = this->nodes[this->get_arc(vertex.incident_edge()->twin()->cell()->source_index()).from_idx].point;
+            const Point &first_point  = this->nodes[this->get_border_arc(vertex.incident_edge()->cell()->source_index()).from_idx].point;
+            const Point &second_point = this->nodes[this->get_border_arc(vertex.incident_edge()->twin()->cell()->source_index()).from_idx].point;
 
             if (vertex_equal_to_point(&vertex, first_point)) {
                 assert(vertex.color() != vertex.incident_edge()->cell()->source_index());
                 assert(vertex.color() != vertex.incident_edge()->twin()->cell()->source_index());
-                vertex.color(this->get_arc(vertex.incident_edge()->cell()->source_index()).from_idx);
+                vertex.color(this->get_border_arc(vertex.incident_edge()->cell()->source_index()).from_idx);
             } else if (vertex_equal_to_point(&vertex, second_point)) {
                 assert(vertex.color() != vertex.incident_edge()->cell()->source_index());
                 assert(vertex.color() != vertex.incident_edge()->twin()->cell()->source_index());
-                vertex.color(this->get_arc(vertex.incident_edge()->twin()->cell()->source_index()).from_idx);
+                vertex.color(this->get_border_arc(vertex.incident_edge()->twin()->cell()->source_index()).from_idx);
             } else if (bbox.contains(vertex_point)) {
                 if (auto [contour_pt, c_dist_sqr] = closest_contour_point.find(vertex_point); contour_pt != nullptr && c_dist_sqr < 3 * SCALED_EPSILON) {
                     vertex.color(this->get_global_index(contour_pt->m_contour_idx, contour_pt->m_point_idx));
@@ -684,6 +708,35 @@ struct MMU_Graph
             }
         }
     }
+
+    void garbage_collect()
+    {
+        std::vector<int> nodes_map(this->nodes.size(), -1);
+        int              nodes_count = 0;
+        size_t           arcs_count  = 0;
+        for (const MMU_Graph::Node &node : this->nodes)
+            if (size_t node_idx = &node - &this->nodes.front(); !node.arc_idxs.empty()) {
+                nodes_map[node_idx] = nodes_count++;
+                arcs_count += node.arc_idxs.size();
+            }
+
+        std::vector<MMU_Graph::Node> new_nodes;
+        std::vector<MMU_Graph::Arc>  new_arcs;
+        new_nodes.reserve(nodes_count);
+        new_arcs.reserve(arcs_count);
+        for (const MMU_Graph::Node &node : this->nodes)
+            if (size_t node_idx = &node - &this->nodes.front(); nodes_map[node_idx] >= 0) {
+                new_nodes.push_back({node.point});
+                for (const size_t &arc_idx : node.arc_idxs) {
+                    const Arc &arc = this->arcs[arc_idx];
+                    new_nodes.back().arc_idxs.emplace_back(new_arcs.size());
+                    new_arcs.push_back({size_t(nodes_map[arc.from_idx]), size_t(nodes_map[arc.to_idx]), arc.color, arc.type});
+                }
+            }
+
+        this->nodes = std::move(new_nodes);
+        this->arcs  = std::move(new_arcs);
+    }
 };
 
 static inline void mark_processed(const voronoi_diagram<double>::const_edge_iterator &edge_iterator)
@@ -825,7 +878,7 @@ static MMU_Graph build_graph(size_t layer_idx, const std::vector<std::vector<Col
             Point              contour_intersection;
 
             if (line_intersection_with_epsilon(contour_line.line, edge_line, &contour_intersection)) {
-                const MMU_Graph::Arc &graph_arc = graph.get_arc(edge_it->cell()->source_index());
+                const MMU_Graph::Arc &graph_arc = graph.get_border_arc(edge_it->cell()->source_index());
                 const size_t          from_idx  = (edge_it->vertex1() != nullptr) ? edge_it->vertex1()->color() : edge_it->vertex0()->color();
                 size_t                to_idx    = ((contour_line.line.a - contour_intersection).cast<double>().squaredNorm() <
                                  (contour_line.line.b - contour_intersection).cast<double>().squaredNorm()) ?
@@ -859,12 +912,12 @@ static MMU_Graph build_graph(size_t layer_idx, const std::vector<std::vector<Col
                 if (edge_it->vertex1()->color() < graph.nodes_count() && !graph.is_vertex_on_contour(edge_it->vertex1())) {
                     Line contour_line_twin = lines_colored[edge_it->twin()->cell()->source_index()].line;
                     if (line_intersection_with_epsilon(contour_line_twin, edge_line, &intersection)) {
-                        const MMU_Graph::Arc &graph_arc = graph.get_arc(edge_it->twin()->cell()->source_index());
+                        const MMU_Graph::Arc &graph_arc = graph.get_border_arc(edge_it->twin()->cell()->source_index());
                         const size_t          to_idx_l  = is_point_closer_to_beginning_of_line(contour_line_twin, intersection) ? graph_arc.from_idx :
                                                                                                                                   graph_arc.to_idx;
                         graph.append_edge(edge_it->vertex1()->color(), to_idx_l);
                     } else if (line_intersection_with_epsilon(contour_line, edge_line, &intersection)) {
-                        const MMU_Graph::Arc &graph_arc = graph.get_arc(edge_it->cell()->source_index());
+                        const MMU_Graph::Arc &graph_arc = graph.get_border_arc(edge_it->cell()->source_index());
                         const size_t to_idx_l = is_point_closer_to_beginning_of_line(contour_line, intersection) ? graph_arc.from_idx : graph_arc.to_idx;
                         graph.append_edge(edge_it->vertex1()->color(), to_idx_l);
                     }
@@ -912,27 +965,25 @@ static MMU_Graph build_graph(size_t layer_idx, const std::vector<std::vector<Col
                     Line second_part(intersection, real_v1);
 
                     if (!has_same_color(contour_line_prev, colored_line)) {
-                        if (points_inside(contour_line_prev.line, contour_line, first_part.b)) {
-                            graph.append_edge(edge_it->vertex0()->color(), graph.get_arc(edge_it->cell()->source_index()).from_idx);
-                        }
-                        if (points_inside(contour_line_prev.line, contour_line, second_part.b)) {
-                            graph.append_edge(edge_it->vertex1()->color(), graph.get_arc(edge_it->cell()->source_index()).from_idx);
-                        }
+                        if (points_inside(contour_line_prev.line, contour_line, first_part.b))
+                            graph.append_edge(edge_it->vertex0()->color(), graph.get_border_arc(edge_it->cell()->source_index()).from_idx);
+
+                        if (points_inside(contour_line_prev.line, contour_line, second_part.b))
+                            graph.append_edge(edge_it->vertex1()->color(), graph.get_border_arc(edge_it->cell()->source_index()).from_idx);
                     }
                 } else {
-                    const size_t int_point_idx = graph.get_arc(edge_it->cell()->source_index()).to_idx;
+                    const size_t int_point_idx = graph.get_border_arc(edge_it->cell()->source_index()).to_idx;
                     const Point  int_point     = graph.nodes[int_point_idx].point;
 
                     const Line first_part(int_point, real_v0);
                     const Line second_part(int_point, real_v1);
 
                     if (!has_same_color(contour_line_next, colored_line)) {
-                        if (points_inside(contour_line, contour_line_next.line, first_part.b)) {
+                        if (points_inside(contour_line, contour_line_next.line, first_part.b))
                             graph.append_edge(edge_it->vertex0()->color(), int_point_idx);
-                        }
-                        if (points_inside(contour_line, contour_line_next.line, second_part.b)) {
+
+                        if (points_inside(contour_line, contour_line_next.line, second_part.b))
                             graph.append_edge(edge_it->vertex1()->color(), int_point_idx);
-                        }
                     }
                 }
             }
@@ -974,13 +1025,15 @@ static inline Polygon to_polygon(const Lines &lines)
 // It iterates through all nodes on the border between two different colors, and from this point,
 // start selection always left most edges for every node to construct CCW polygons.
 // Assumes that graph is planar (without self-intersection edges)
-static std::vector<std::pair<Polygon, size_t>> extract_colored_segments(MMU_Graph &graph)
+static std::vector<std::pair<Polygon, size_t>> extract_colored_segments(const MMU_Graph &graph)
 {
+    std::vector<bool> used_arcs(graph.arcs.size(), false);
     // When there is no next arc, then is returned original_arc or edge with is marked as used
-    auto get_next = [&graph](const Line &process_line, MMU_Graph::Arc &original_arc) -> MMU_Graph::Arc & {
-        std::vector<std::pair<MMU_Graph::Arc *, double>> sorted_arcs;
-        for (MMU_Graph::Arc &arc : graph.nodes[original_arc.to_idx].neighbours) {
-            if (graph.nodes[arc.to_idx].point == process_line.a || arc.used)
+    auto get_next = [&graph, &used_arcs](const Line &process_line, const MMU_Graph::Arc &original_arc) -> const MMU_Graph::Arc & {
+        std::vector<std::pair<const MMU_Graph::Arc *, double>> sorted_arcs;
+        for (const size_t &arc_idx : graph.nodes[original_arc.to_idx].arc_idxs) {
+            const MMU_Graph::Arc &arc = graph.arcs[arc_idx];
+            if (graph.nodes[arc.to_idx].point == process_line.a || used_arcs[arc_idx])
                 continue;
 
             assert(original_arc.to_idx == arc.from_idx);
@@ -995,11 +1048,11 @@ static std::vector<std::pair<Polygon, size_t>> extract_colored_segments(MMU_Grap
         }
 
         std::sort(sorted_arcs.begin(), sorted_arcs.end(),
-                  [](std::pair<MMU_Graph::Arc *, double> &l, std::pair<MMU_Graph::Arc *, double> &r) -> bool { return l.second < r.second; });
+                  [](std::pair<const MMU_Graph::Arc *, double> &l, std::pair<const MMU_Graph::Arc *, double> &r) -> bool { return l.second < r.second; });
 
         // Try to return left most edge witch is unused
         for (auto &sorted_arc : sorted_arcs)
-            if (!sorted_arc.first->used)
+            if (size_t arc_idx = sorted_arc.first - &graph.arcs.front(); !used_arcs[arc_idx])
                 return *sorted_arc.first;
 
         if (sorted_arcs.empty())
@@ -1008,35 +1061,39 @@ static std::vector<std::pair<Polygon, size_t>> extract_colored_segments(MMU_Grap
         return *(sorted_arcs.front().first);
     };
 
+    auto all_arc_used = [&used_arcs](const MMU_Graph::Node &node) -> bool {
+        return std::all_of(node.arc_idxs.cbegin(), node.arc_idxs.cend(), [&used_arcs](const size_t &arc_idx) -> bool { return used_arcs[arc_idx]; });
+    };
+
     std::vector<std::pair<Polygon, size_t>> polygons_segments;
-    for (MMU_Graph::Node &node : graph.nodes)
-        for (MMU_Graph::Arc &arc : node.neighbours)
-            arc.used = false;
-
     for (size_t node_idx = 0; node_idx < graph.all_border_points; ++node_idx) {
-        MMU_Graph::Node &node = graph.nodes[node_idx];
+        const MMU_Graph::Node &node = graph.nodes[node_idx];
+
+        for (const size_t &arc_idx : node.arc_idxs) {
+            const MMU_Graph::Arc &arc = graph.arcs[arc_idx];
+            if (arc.type == MMU_Graph::ARC_TYPE::NON_BORDER || used_arcs[arc_idx])continue;
 
-        for (MMU_Graph::Arc &arc : node.neighbours) {
-            if (arc.type == MMU_Graph::ARC_TYPE::NON_BORDER || arc.used) continue;
 
             Line process_line(node.point, graph.nodes[arc.to_idx].point);
-            arc.used = true;
+            used_arcs[arc_idx] = true;
 
             Lines face_lines;
             face_lines.emplace_back(process_line);
             Point start_p = process_line.a;
 
-            Line            p_vec = process_line;
-            MMU_Graph::Arc *p_arc = &arc;
+            Line                  p_vec = process_line;
+            const MMU_Graph::Arc *p_arc = &arc;
             do {
-                MMU_Graph::Arc &next = get_next(p_vec, *p_arc);
-                face_lines.emplace_back(Line(graph.nodes[next.from_idx].point, graph.nodes[next.to_idx].point));
-                if (next.used) break;
+                const MMU_Graph::Arc &next         = get_next(p_vec, *p_arc);
+                size_t                next_arc_idx = &next - &graph.arcs.front();
+                face_lines.emplace_back(graph.nodes[next.from_idx].point, graph.nodes[next.to_idx].point);
+                if (used_arcs[next_arc_idx])
+                    break;
 
-                next.used = true;
-                p_vec     = Line(graph.nodes[next.from_idx].point, graph.nodes[next.to_idx].point);
-                p_arc     = &next;
-            } while (graph.nodes[p_arc->to_idx].point != start_p);
+                used_arcs[next_arc_idx] = true;
+                p_vec                   = Line(graph.nodes[next.from_idx].point, graph.nodes[next.to_idx].point);
+                p_arc                   = &next;
+            } while (graph.nodes[p_arc->to_idx].point != start_p || !all_arc_used(graph.nodes[p_arc->to_idx]));
 
             Polygon poly = to_polygon(face_lines);
             if (poly.is_counter_clockwise() && poly.is_valid())
@@ -1049,20 +1106,19 @@ static std::vector<std::pair<Polygon, size_t>> extract_colored_segments(MMU_Grap
 // Used in remove_multiple_edges_in_vertices()
 // Returns length of edge with is connected to contour. To this length is include other edges with follows it if they are almost straight (with the
 // tolerance of 15) And also if node between two subsequent edges is connected only to these two edges.
-static inline double compute_edge_length(MMU_Graph &graph, size_t start_idx, MMU_Graph::Arc &start_edge)
+static inline double compute_edge_length(const MMU_Graph &graph, const size_t start_idx, const size_t &start_arc_idx)
 {
-    for (MMU_Graph::Node &node : graph.nodes)
-        for (MMU_Graph::Arc &arc : node.neighbours)
-            arc.used = false;
+    assert(start_arc_idx < graph.arcs.size());
+    std::vector<bool> used_arcs(graph.arcs.size(), false);
 
-    start_edge.used                   = true;
-    MMU_Graph::Arc *arc               = &start_edge;
-    size_t          idx               = start_idx;
-    double          line_total_length = Line(graph.nodes[idx].point, graph.nodes[arc->to_idx].point).length();
-    while (graph.nodes[arc->to_idx].neighbours.size() == 2) {
+    used_arcs[start_arc_idx]                = true;
+    const MMU_Graph::Arc *arc               = &graph.arcs[start_arc_idx];
+    size_t                idx               = start_idx;
+    double                line_total_length = (graph.nodes[arc->to_idx].point - graph.nodes[idx].point).cast<double>().norm();;
+    while (graph.nodes[arc->to_idx].arc_idxs.size() == 2) {
         bool found = false;
-        for (MMU_Graph::Arc &arc_n : graph.nodes[arc->to_idx].neighbours) {
-            if (arc_n.type == MMU_Graph::ARC_TYPE::NON_BORDER && !arc_n.used && arc_n.to_idx != idx) {
+        for (const size_t &arc_idx : graph.nodes[arc->to_idx].arc_idxs) {
+            if (const MMU_Graph::Arc &arc_n = graph.arcs[arc_idx]; arc_n.type == MMU_Graph::ARC_TYPE::NON_BORDER && !used_arcs[arc_idx] && arc_n.to_idx != idx) {
                 Line first_line(graph.nodes[idx].point, graph.nodes[arc->to_idx].point);
                 Line second_line(graph.nodes[arc->to_idx].point, graph.nodes[arc_n.to_idx].point);
 
@@ -1080,8 +1136,8 @@ static inline double compute_edge_length(MMU_Graph &graph, size_t start_idx, MMU
                 idx = arc->to_idx;
                 arc = &arc_n;
 
-                line_total_length += Line(graph.nodes[idx].point, graph.nodes[arc->to_idx].point).length();
-                arc_n.used = true;
+                line_total_length += (graph.nodes[arc->to_idx].point - graph.nodes[idx].point).cast<double>().norm();
+                used_arcs[arc_idx] = true;
                 found      = true;
                 break;
             }
@@ -1104,11 +1160,12 @@ static void remove_multiple_edges_in_vertices(MMU_Graph &graph, const std::vecto
             size_t second_idx = graph.get_global_index(poly_idx, (colored_segment.second + 1) % graph.polygon_sizes[poly_idx]);
             Line   seg_line(graph.nodes[first_idx].point, graph.nodes[second_idx].point);
 
-            if (graph.nodes[first_idx].neighbours.size() >= 3) {
+            if (graph.nodes[first_idx].arc_idxs.size() >= 3) {
                 std::vector<std::pair<MMU_Graph::Arc *, double>> arc_to_check;
-                for (MMU_Graph::Arc &n_arc : graph.nodes[first_idx].neighbours) {
+                for (const size_t &arc_idx : graph.nodes[first_idx].arc_idxs) {
+                    MMU_Graph::Arc &n_arc = graph.arcs[arc_idx];
                     if (n_arc.type == MMU_Graph::ARC_TYPE::NON_BORDER) {
-                        double total_len = compute_edge_length(graph, first_idx, n_arc);
+                        double total_len = compute_edge_length(graph, first_idx, arc_idx);
                         arc_to_check.emplace_back(&n_arc, total_len);
                     }
                 }
@@ -1478,18 +1535,18 @@ static inline std::vector<std::vector<ExPolygons>> mmu_segmentation_top_and_bott
         LayerColorStat out;
         const Layer &layer = *layers[layer_idx];
         for (const LayerRegion *region : layer.regions())
-            if (const PrintRegionConfig &config = region->region().config(); 
+            if (const PrintRegionConfig &config = region->region().config();
                 // color_idx == 0 means "don't know" extruder aka the underlying extruder.
                 // As this region may split existing regions, we collect statistics over all regions for color_idx == 0.
                 color_idx == 0 || config.perimeter_extruder == int(color_idx)) {
-                out.extrusion_width     = std::max<float>(out.extrusion_width, config.perimeter_extrusion_width);
-                out.top_solid_layers    = std::max<float>(out.top_solid_layers, config.top_solid_layers);
-                out.bottom_solid_layers = std::max<float>(out.bottom_solid_layers, config.bottom_solid_layers);
-                out.small_region_threshold = config.gap_fill_enabled.value && config.gap_fill_speed.value > 0 ? 
-                    // Gap fill enabled. Enable a single line of 1/2 extrusion width.
-                    0.5 * config.perimeter_extrusion_width :
-                    // Gap fill disabled. Enable two lines slightly overlapping.
-                    config.perimeter_extrusion_width + 0.7f * Flow::rounded_rectangle_extrusion_spacing(config.perimeter_extrusion_width, layer.height);
+                out.extrusion_width     = std::max<float>(out.extrusion_width, float(config.perimeter_extrusion_width));
+                out.top_solid_layers    = std::max<int>(out.top_solid_layers, config.top_solid_layers);
+                out.bottom_solid_layers = std::max<int>(out.bottom_solid_layers, config.bottom_solid_layers);
+                out.small_region_threshold = config.gap_fill_enabled.value && config.gap_fill_speed.value > 0 ?
+                                             // Gap fill enabled. Enable a single line of 1/2 extrusion width.
+                                             0.5f * float(config.perimeter_extrusion_width) :
+                                             // Gap fill disabled. Enable two lines slightly overlapping.
+                                             float(config.perimeter_extrusion_width) + 0.7f * Flow::rounded_rectangle_extrusion_spacing(float(config.perimeter_extrusion_width), float(layer.height));
                 out.small_region_threshold = scaled<float>(out.small_region_threshold * 0.5f);
                 ++ out.num_regions;
             }
@@ -1603,14 +1660,70 @@ static std::vector<std::vector<std::pair<ExPolygon, size_t>>> merge_segmented_la
     return segmented_regions_merged;
 }
 
+#ifdef MMU_SEGMENTATION_DEBUG_REGIONS
+static void export_regions_to_svg(const std::string &path, const std::vector<std::pair<ExPolygon, size_t>> &regions, const ExPolygons &lslices)
+{
+    const std::vector<std::string> colors       = {"blue", "cyan", "red", "orange", "magenta", "pink", "purple", "yellow"};
+    coordf_t                       stroke_width = scale_(0.05);
+    BoundingBox                    bbox         = get_extents(lslices);
+    bbox.offset(scale_(1.));
+    ::Slic3r::SVG svg(path.c_str(), bbox);
+
+    svg.draw_outline(lslices, "green", "lime", stroke_width);
+    for (const std::pair<ExPolygon, size_t> &region : regions) {
+        int region_color = region.second;
+        if (region_color >= 0 && region_color < int(colors.size()))
+            svg.draw(region.first, colors[region_color]);
+        else
+            svg.draw(region.first, "black");
+    }
+}
+#endif // MMU_SEGMENTATION_DEBUG_REGIONS
+
+#ifdef MMU_SEGMENTATION_DEBUG_GRAPH
+static void export_graph_to_svg(const std::string &path, const MMU_Graph &graph, const ExPolygons &lslices)
+{
+    const std::vector<std::string> colors       = {"blue", "cyan", "red", "orange", "magenta", "pink", "purple", "green", "yellow"};
+    coordf_t                       stroke_width = scale_(0.05);
+    BoundingBox                    bbox         = get_extents(lslices);
+    bbox.offset(scale_(1.));
+    ::Slic3r::SVG svg(path.c_str(), bbox);
+    for (const MMU_Graph::Node &node : graph.nodes)
+        for (const size_t &arc_idx : node.arc_idxs) {
+            const MMU_Graph::Arc &arc = graph.arcs[arc_idx];
+            Line arc_line(node.point, graph.nodes[arc.to_idx].point);
+            if (arc.type == MMU_Graph::ARC_TYPE::BORDER && arc.color >= 0 && arc.color < int(colors.size()))
+                svg.draw(arc_line, colors[arc.color], stroke_width);
+            else
+                svg.draw(arc_line, "black", stroke_width);
+        }
+}
+#endif // MMU_SEGMENTATION_DEBUG_GRAPH
+
+#ifdef MMU_SEGMENTATION_DEBUG_INPUT
+void export_processed_input_expolygons_to_svg(const std::string &path, const LayerRegionPtrs &regions, const ExPolygons &processed_input_expolygons)
+{
+    coordf_t    stroke_width = scale_(0.05);
+    BoundingBox bbox         = get_extents(regions);
+    bbox.merge(get_extents(processed_input_expolygons));
+    bbox.offset(scale_(1.));
+    ::Slic3r::SVG svg(path.c_str(), bbox);
+
+    for (LayerRegion *region : regions)
+        svg.draw_outline(region->slices.surfaces, "blue", "cyan", stroke_width);
+
+    svg.draw_outline(processed_input_expolygons, "red", "pink", stroke_width);
+}
+#endif // MMU_SEGMENTATION_DEBUG_INPUT
+
 std::vector<std::vector<std::pair<ExPolygon, size_t>>> multi_material_segmentation_by_painting(const PrintObject &print_object, const std::function<void()> &throw_on_cancel_callback)
 {
     std::vector<std::vector<std::pair<ExPolygon, size_t>>> segmented_regions(print_object.layers().size());
     std::vector<std::vector<PaintedLine>>                  painted_lines(print_object.layers().size());
+    std::array<std::mutex, 64>                             painted_lines_mutex;
     std::vector<EdgeGrid::Grid>                            edge_grids(print_object.layers().size());
     const ConstLayerPtrsAdaptor                            layers = print_object.layers();
     std::vector<ExPolygons>                                input_expolygons(layers.size());
-    std::vector<Polygons>                                  input_polygons(layers.size());
 
     throw_on_cancel_callback();
 
@@ -1636,86 +1749,99 @@ std::vector<std::vector<std::pair<ExPolygon, size_t>>> multi_material_segmentati
             // This consequently leads to issues with the extraction of colored segments by function extract_colored_segments.
             // Calling expolygons_simplify fixed these issues.
             input_expolygons[layer_idx] = smooth_outward(expolygons_simplify(offset_ex(ex_polygons, -10.f * float(SCALED_EPSILON)), 5 * SCALED_EPSILON), 10 * coord_t(SCALED_EPSILON));
-            input_polygons[layer_idx]   = to_polygons(input_expolygons[layer_idx]);
+
+#ifdef MMU_SEGMENTATION_DEBUG_INPUT
+            {
+                static int iRun = 0;
+                export_processed_input_expolygons_to_svg(debug_out_path("mm-input-%d-%d.svg", layer_idx, iRun++), layers[layer_idx]->regions(), input_expolygons[layer_idx]);
+            }
+#endif // MMU_SEGMENTATION_DEBUG_INPUT
         }
     }); // end of parallel_for
     BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - slices preparation in parallel - end";
 
     for (size_t layer_idx = 0; layer_idx < layers.size(); ++layer_idx) {
         throw_on_cancel_callback();
-        BoundingBox  bbox(get_extents(input_polygons[layer_idx]));
+        BoundingBox bbox(get_extents(layers[layer_idx]->regions()));
+        bbox.merge(get_extents(input_expolygons[layer_idx]));
         // Projected triangles may slightly exceed the input polygons.
         bbox.offset(20 * SCALED_EPSILON);
         edge_grids[layer_idx].set_bbox(bbox);
-        edge_grids[layer_idx].create(input_polygons[layer_idx], coord_t(scale_(10.)));
+        edge_grids[layer_idx].create(input_expolygons[layer_idx], coord_t(scale_(10.)));
     }
 
     BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - projection of painted triangles - begin";
     for (const ModelVolume *mv : print_object.model_object()->volumes) {
         const size_t num_extruders = print_object.print()->config().nozzle_diameter.size() + 1;
-        for (size_t extruder_idx = 1; extruder_idx < num_extruders; ++extruder_idx) {
-            throw_on_cancel_callback();
-            const indexed_triangle_set custom_facets = mv->mmu_segmentation_facets.get_facets(*mv, EnforcerBlockerType(extruder_idx));
-            if (!mv->is_model_part() || custom_facets.indices.empty())
-                continue;
+        tbb::parallel_for(tbb::blocked_range<size_t>(1, num_extruders), [&mv, &print_object, &edge_grids, &painted_lines, &painted_lines_mutex, &input_expolygons, &throw_on_cancel_callback](const tbb::blocked_range<size_t> &range) {
+            for (size_t extruder_idx = range.begin(); extruder_idx < range.end(); ++extruder_idx) {
+                throw_on_cancel_callback();
+                const indexed_triangle_set custom_facets = mv->mmu_segmentation_facets.get_facets(*mv, EnforcerBlockerType(extruder_idx));
+                if (!mv->is_model_part() || custom_facets.indices.empty())
+                    continue;
 
-            const Transform3f tr = print_object.trafo().cast<float>() * mv->get_matrix().cast<float>();
-            for (size_t facet_idx = 0; facet_idx < custom_facets.indices.size(); ++facet_idx) {
-                float min_z = std::numeric_limits<float>::max();
-                float max_z = std::numeric_limits<float>::lowest();
+                const Transform3f tr = print_object.trafo().cast<float>() * mv->get_matrix().cast<float>();
+                tbb::parallel_for(tbb::blocked_range<size_t>(0, custom_facets.indices.size()), [&tr, &custom_facets, &print_object, &edge_grids, &input_expolygons, &painted_lines, &painted_lines_mutex, &extruder_idx](const tbb::blocked_range<size_t> &range) {
+                    for (size_t facet_idx = range.begin(); facet_idx < range.end(); ++facet_idx) {
+                        float min_z = std::numeric_limits<float>::max();
+                        float max_z = std::numeric_limits<float>::lowest();
 
-                std::array<Vec3f, 3> facet;
-                for (int p_idx = 0; p_idx < 3; ++p_idx) {
-                    facet[p_idx] = tr * custom_facets.vertices[custom_facets.indices[facet_idx](p_idx)];
-                    max_z        = std::max(max_z, facet[p_idx].z());
-                    min_z        = std::min(min_z, facet[p_idx].z());
-                }
+                        std::array<Vec3f, 3> facet;
+                        for (int p_idx = 0; p_idx < 3; ++p_idx) {
+                            facet[p_idx] = tr * custom_facets.vertices[custom_facets.indices[facet_idx](p_idx)];
+                            max_z        = std::max(max_z, facet[p_idx].z());
+                            min_z        = std::min(min_z, facet[p_idx].z());
+                        }
 
-                // Sort the vertices by z-axis for simplification of projected_facet on slices
-                std::sort(facet.begin(), facet.end(), [](const Vec3f &p1, const Vec3f &p2) { return p1.z() < p2.z(); });
+                        // Sort the vertices by z-axis for simplification of projected_facet on slices
+                        std::sort(facet.begin(), facet.end(), [](const Vec3f &p1, const Vec3f &p2) { return p1.z() < p2.z(); });
 
-                // Find lowest slice not below the triangle.
-                auto first_layer = std::upper_bound(print_object.layers().begin(), print_object.layers().end(), float(min_z - EPSILON),
-                                                    [](float z, const Layer *l1) { return z < l1->slice_z; });
-                auto last_layer  = std::upper_bound(print_object.layers().begin(), print_object.layers().end(), float(max_z + EPSILON),
-                                                   [](float z, const Layer *l1) { return z < l1->slice_z; });
-                --last_layer;
+                        // Find lowest slice not below the triangle.
+                        auto first_layer = std::upper_bound(print_object.layers().begin(), print_object.layers().end(), float(min_z - EPSILON),
+                                                            [](float z, const Layer *l1) { return z < l1->slice_z; });
+                        auto last_layer  = std::upper_bound(print_object.layers().begin(), print_object.layers().end(), float(max_z + EPSILON),
+                                                           [](float z, const Layer *l1) { return z < l1->slice_z; });
+                        --last_layer;
 
-                for (auto layer_it = first_layer; layer_it != (last_layer + 1); ++layer_it) {
-                    const Layer *layer     = *layer_it;
-                    size_t       layer_idx = layer_it - print_object.layers().begin();
-                    if (facet[0].z() > layer->slice_z || layer->slice_z > facet[2].z())
-                        continue;
+                        for (auto layer_it = first_layer; layer_it != (last_layer + 1); ++layer_it) {
+                            const Layer *layer     = *layer_it;
+                            size_t       layer_idx = layer_it - print_object.layers().begin();
+                            if (input_expolygons[layer_idx].empty() || facet[0].z() > layer->slice_z || layer->slice_z > facet[2].z())
+                                continue;
 
-                    // https://kandepet.com/3d-printing-slicing-3d-objects/
-                    float t            = (float(layer->slice_z) - facet[0].z()) / (facet[2].z() - facet[0].z());
-                    Vec3f line_start_f = facet[0] + t * (facet[2] - facet[0]);
-                    Vec3f line_end_f;
+                            // https://kandepet.com/3d-printing-slicing-3d-objects/
+                            float t            = (float(layer->slice_z) - facet[0].z()) / (facet[2].z() - facet[0].z());
+                            Vec3f line_start_f = facet[0] + t * (facet[2] - facet[0]);
+                            Vec3f line_end_f;
 
-                    if (facet[1].z() > layer->slice_z) {
-                        // [P0, P2] a [P0, P1]
-                        float t1   = (float(layer->slice_z) - facet[0].z()) / (facet[1].z() - facet[0].z());
-                        line_end_f = facet[0] + t1 * (facet[1] - facet[0]);
-                    } else {
-                        // [P0, P2] a [P1, P2]
-                        float t2   = (float(layer->slice_z) - facet[1].z()) / (facet[2].z() - facet[1].z());
-                        line_end_f = facet[1] + t2 * (facet[2] - facet[1]);
+                            if (facet[1].z() > layer->slice_z) {
+                                // [P0, P2] and [P0, P1]
+                                float t1   = (float(layer->slice_z) - facet[0].z()) / (facet[1].z() - facet[0].z());
+                                line_end_f = facet[0] + t1 * (facet[1] - facet[0]);
+                            } else {
+                                // [P0, P2] and [P1, P2]
+                                float t2   = (float(layer->slice_z) - facet[1].z()) / (facet[2].z() - facet[1].z());
+                                line_end_f = facet[1] + t2 * (facet[2] - facet[1]);
+                            }
+
+                            Point line_start(scale_(line_start_f.x()), scale_(line_start_f.y()));
+                            Point line_end(scale_(line_end_f.x()), scale_(line_end_f.y()));
+                            line_start -= print_object.center_offset();
+                            line_end   -= print_object.center_offset();
+
+                            size_t mutex_idx = layer_idx & 0x3F;
+                            assert(mutex_idx < painted_lines_mutex.size());
+
+                            PaintedLineVisitor visitor(edge_grids[layer_idx], painted_lines[layer_idx], painted_lines_mutex[mutex_idx], 16);
+                            visitor.line_to_test.a = line_start;
+                            visitor.line_to_test.b = line_end;
+                            visitor.color          = int(extruder_idx);
+                            edge_grids[layer_idx].visit_cells_intersecting_line(line_start, line_end, visitor);
+                        }
                     }
-
-                    Point line_start(scale_(line_start_f.x()), scale_(line_start_f.y()));
-                    Point line_end(scale_(line_end_f.x()), scale_(line_end_f.y()));
-                    line_start -= print_object.center_offset();
-                    line_end   -= print_object.center_offset();
-
-                    PaintedLineVisitor       visitor(edge_grids[layer_idx], painted_lines[layer_idx], 16);
-                    visitor.reset();
-                    visitor.line_to_test.a = line_start;
-                    visitor.line_to_test.b = line_end;
-                    visitor.color          = int(extruder_idx);
-                    edge_grids[layer_idx].visit_cells_intersecting_line(line_start, line_end, visitor);
-                }
+                });
             }
-        }
+        });
     }
     BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - projection of painted triangles - end";
     BOOST_LOG_TRIVIAL(debug) << "MMU segmentation - painted layers count: "
@@ -1725,8 +1851,8 @@ std::vector<std::vector<std::pair<ExPolygon, size_t>>> multi_material_segmentati
     tbb::parallel_for(tbb::blocked_range<size_t>(0, print_object.layers().size()), [&](const tbb::blocked_range<size_t> &range) {
         for (size_t layer_idx = range.begin(); layer_idx < range.end(); ++layer_idx) {
             throw_on_cancel_callback();
-            auto comp = [&input_polygons, layer_idx](const PaintedLine &first, const PaintedLine &second) {
-                Point first_start_p = input_polygons[layer_idx][first.contour_idx][first.line_idx];
+            auto comp = [&edge_grids, layer_idx](const PaintedLine &first, const PaintedLine &second) {
+                Point first_start_p = edge_grids[layer_idx].contours()[first.contour_idx].segment_start(first.line_idx);
                 return first.contour_idx < second.contour_idx ||
                        (first.contour_idx == second.contour_idx &&
                         (first.line_idx < second.line_idx ||
@@ -1740,13 +1866,28 @@ std::vector<std::vector<std::pair<ExPolygon, size_t>>> multi_material_segmentati
             std::vector<PaintedLine> &painted_lines_single = painted_lines[layer_idx];
 
             if (!painted_lines_single.empty()) {
-                std::vector<std::vector<ColoredLine>> color_poly = colorize_polygons(input_polygons[layer_idx], painted_lines_single);
+                std::vector<std::vector<ColoredLine>> color_poly = colorize_polygons(edge_grids[layer_idx].contours(), painted_lines_single);
                 MMU_Graph                             graph      = build_graph(layer_idx, color_poly);
                 remove_multiple_edges_in_vertices(graph, color_poly);
                 graph.remove_nodes_with_one_arc();
+
+#ifdef MMU_SEGMENTATION_DEBUG_GRAPH
+                {
+                    static int iRun = 0;
+                    export_graph_to_svg(debug_out_path("mm-graph-final-%d-%d.svg", layer_idx, iRun++), graph, input_expolygons[layer_idx]);
+                }
+#endif // MMU_SEGMENTATION_DEBUG_GRAPH
+
                 std::vector<std::pair<Polygon, size_t>> segmentation = extract_colored_segments(graph);
                 for (std::pair<Polygon, size_t> &region : segmentation)
                     segmented_regions[layer_idx].emplace_back(std::move(region));
+
+#ifdef MMU_SEGMENTATION_DEBUG_REGIONS
+                {
+                    static int iRun = 0;
+                    export_regions_to_svg(debug_out_path("mm-regions-sides-%d-%d.svg", layer_idx, iRun++), segmented_regions[layer_idx], input_expolygons[layer_idx]);
+                }
+#endif // MMU_SEGMENTATION_DEBUG_REGIONS
             }
         }
     }); // end of parallel_for
@@ -1765,6 +1906,14 @@ std::vector<std::vector<std::pair<ExPolygon, size_t>>> multi_material_segmentati
     std::vector<std::vector<std::pair<ExPolygon, size_t>>> segmented_regions_merged = merge_segmented_layers(segmented_regions, std::move(top_and_bottom_layers), throw_on_cancel_callback);
     throw_on_cancel_callback();
 
+#ifdef MMU_SEGMENTATION_DEBUG_REGIONS
+    {
+        static int iRun = 0;
+        for (size_t layer_idx = 0; layer_idx < print_object.layers().size(); ++layer_idx)
+            export_regions_to_svg(debug_out_path("mm-regions-merged-%d-%d.svg", layer_idx, iRun++), segmented_regions_merged[layer_idx], input_expolygons[layer_idx]);
+    }
+#endif // MMU_SEGMENTATION_DEBUG_REGIONS
+
     return segmented_regions_merged;
 }
 
diff --git a/src/libslic3r/PrintConfig.cpp b/src/libslic3r/PrintConfig.cpp
index 3c5ff3859..02ce4fa4b 100644
--- a/src/libslic3r/PrintConfig.cpp
+++ b/src/libslic3r/PrintConfig.cpp
@@ -1457,6 +1457,7 @@ void PrintConfigDef::init_fff_params()
     def->tooltip = L("Maximum width of a segmented region. Zero disables this feature.");
     def->sidetext = L("mm (zero to disable)");
     def->min = 0;
+    def->category = L("Advanced");
     def->mode = comExpert;
     def->set_default_value(new ConfigOptionFloat(0.f));
 
diff --git a/src/libslic3r/Technologies.hpp b/src/libslic3r/Technologies.hpp
index 339082b52..45217c959 100644
--- a/src/libslic3r/Technologies.hpp
+++ b/src/libslic3r/Technologies.hpp
@@ -10,8 +10,6 @@
 #define ENABLE_SELECTION_DEBUG_OUTPUT 0
 // Renders a small sphere in the center of the bounding box of the current selection when no gizmo is active
 #define ENABLE_RENDER_SELECTION_CENTER 0
-// Shows an imgui dialog with render related data
-#define ENABLE_RENDER_STATISTICS 0
 // Shows an imgui dialog with camera related data
 #define ENABLE_CAMERA_STATISTICS 0
 // Render the picking pass instead of the main scene (use [T] key to toggle between regular rendering and picking pass only rendering)
diff --git a/src/libslic3r/TriangleSelector.cpp b/src/libslic3r/TriangleSelector.cpp
index a349fe35a..ad823c55d 100644
--- a/src/libslic3r/TriangleSelector.cpp
+++ b/src/libslic3r/TriangleSelector.cpp
@@ -228,7 +228,7 @@ void TriangleSelector::seed_fill_select_triangles(const Vec3f &hit, int facet_st
     }
 }
 
-void TriangleSelector::precompute_all_level_neighbors_recursive(const int facet_idx, const Vec3i &neighbors, const Vec3i &neighbors_propagated, std::vector<Vec3i> &neighbors_out) const
+void TriangleSelector::precompute_all_neighbors_recursive(const int facet_idx, const Vec3i &neighbors, const Vec3i &neighbors_propagated, std::vector<Vec3i> &neighbors_out, std::vector<Vec3i> &neighbors_propagated_out) const
 {
     assert(facet_idx < int(m_triangles.size()));
 
@@ -236,7 +236,8 @@ void TriangleSelector::precompute_all_level_neighbors_recursive(const int facet_
     if (!tr->valid())
         return;
 
-    neighbors_out[facet_idx] = neighbors_propagated;
+    neighbors_out[facet_idx]            = neighbors;
+    neighbors_propagated_out[facet_idx] = neighbors_propagated;
     if (tr->is_split()) {
         assert(this->verify_triangle_neighbors(*tr, neighbors));
 
@@ -247,67 +248,51 @@ void TriangleSelector::precompute_all_level_neighbors_recursive(const int facet_
                 assert(tr->children[i] < int(m_triangles.size()));
                 // Recursion, deep first search over the children of this triangle.
                 // All children of this triangle were created by splitting a single source triangle of the original mesh.
-                this->precompute_all_level_neighbors_recursive(tr->children[i], this->child_neighbors(*tr, neighbors, i), this->child_neighbors_propagated(*tr, neighbors_propagated, i), neighbors_out);
+                this->precompute_all_neighbors_recursive(tr->children[i], this->child_neighbors(*tr, neighbors, i),
+                                                         this->child_neighbors_propagated(*tr, neighbors_propagated, i), neighbors_out,
+                                                         neighbors_propagated_out);
             }
         }
     }
 }
 
-std::vector<Vec3i> TriangleSelector::precompute_all_level_neighbors() const
+std::pair<std::vector<Vec3i>, std::vector<Vec3i>> TriangleSelector::precompute_all_neighbors() const
 {
     std::vector<Vec3i> neighbors(m_triangles.size(), Vec3i(-1, -1, -1));
+    std::vector<Vec3i> neighbors_propagated(m_triangles.size(), Vec3i(-1, -1, -1));
     for (int facet_idx = 0; facet_idx < this->m_orig_size_indices; ++facet_idx) {
-        neighbors[facet_idx] = root_neighbors(*m_mesh, facet_idx);
+        neighbors[facet_idx]            = root_neighbors(*m_mesh, facet_idx);
+        neighbors_propagated[facet_idx] = neighbors[facet_idx];
         assert(this->verify_triangle_neighbors(m_triangles[facet_idx], neighbors[facet_idx]));
         if (m_triangles[facet_idx].is_split())
-            this->precompute_all_level_neighbors_recursive(facet_idx, neighbors[facet_idx], neighbors[facet_idx], neighbors);
+            this->precompute_all_neighbors_recursive(facet_idx, neighbors[facet_idx], neighbors_propagated[facet_idx], neighbors, neighbors_propagated);
     }
-    return neighbors;
+    return std::make_pair(std::move(neighbors), std::move(neighbors_propagated));
 }
 
-bool TriangleSelector::are_triangles_touching(const int first_facet_idx, const int second_facet_idx) const
+// It appends all triangles that are touching the edge (vertexi, vertexj) of the triangle.
+// It doesn't append the triangles that are touching the triangle only by part of the edge that means the triangles are from lower depth.
+void TriangleSelector::append_touching_subtriangles(int itriangle, int vertexi, int vertexj, std::vector<int> &touching_subtriangles_out) const
 {
-    std::array<Linef3, 3> sides_facet = {Linef3(m_vertices[m_triangles[first_facet_idx].verts_idxs[0]].v.cast<double>(), m_vertices[m_triangles[first_facet_idx].verts_idxs[1]].v.cast<double>()),
-                                         Linef3(m_vertices[m_triangles[first_facet_idx].verts_idxs[1]].v.cast<double>(), m_vertices[m_triangles[first_facet_idx].verts_idxs[2]].v.cast<double>()),
-                                         Linef3(m_vertices[m_triangles[first_facet_idx].verts_idxs[2]].v.cast<double>(), m_vertices[m_triangles[first_facet_idx].verts_idxs[0]].v.cast<double>())};
+    if (itriangle == -1)
+        return;
 
-    const Vec3d           p0          = m_vertices[m_triangles[second_facet_idx].verts_idxs[0]].v.cast<double>();
-    const Vec3d           p1          = m_vertices[m_triangles[second_facet_idx].verts_idxs[1]].v.cast<double>();
-    const Vec3d           p2          = m_vertices[m_triangles[second_facet_idx].verts_idxs[2]].v.cast<double>();
+    auto process_subtriangle = [this, &itriangle, &vertexi, &vertexj, &touching_subtriangles_out](const int subtriangle_idx) -> void {
+        assert(subtriangle_idx == -1);
+        if (!m_triangles[subtriangle_idx].is_split())
+            touching_subtriangles_out.emplace_back(subtriangle_idx);
+        else if (int midpoint = this->triangle_midpoint(itriangle, vertexi, vertexj); midpoint != -1)
+            append_touching_subtriangles(subtriangle_idx, vertexi, midpoint, touching_subtriangles_out);
+        else
+            append_touching_subtriangles(subtriangle_idx, vertexi, vertexj, touching_subtriangles_out);
+    };
 
-    for (size_t idx = 0; idx < 3; ++idx)
-        if (line_alg::distance_to_squared(sides_facet[idx], p0) <= EPSILON && (line_alg::distance_to_squared(sides_facet[idx], p1) <= EPSILON || line_alg::distance_to_squared(sides_facet[idx], p2) <= EPSILON))
-            return true;
-        else if (line_alg::distance_to_squared(sides_facet[idx], p1) <= EPSILON && line_alg::distance_to_squared(sides_facet[idx], p2) <= EPSILON)
-            return true;
+    std::pair<int, int> touching = this->triangle_subtriangles(itriangle, vertexi, vertexj);
+    if (touching.first != -1)
+        process_subtriangle(touching.first);
 
-    return false;
-}
-
-std::vector<int> TriangleSelector::neighboring_triangles(const int first_facet_idx, const int second_facet_idx, EnforcerBlockerType second_facet_state) const
-{
-    assert(first_facet_idx < int(m_triangles.size()));
-
-    const Triangle *tr = &m_triangles[first_facet_idx];
-    if (!tr->valid())
-        return {};
-
-    if (!tr->is_split() && tr->get_state() == second_facet_state && (are_triangles_touching(second_facet_idx, first_facet_idx) || are_triangles_touching(first_facet_idx, second_facet_idx)))
-        return {first_facet_idx};
-
-    std::vector<int> neighbor_facets_out;
-    int              num_of_children = tr->number_of_split_sides() + 1;
-    if (num_of_children != 1) {
-        for (int i = 0; i < num_of_children; ++i) {
-            assert(i < int(tr->children.size()));
-            assert(tr->children[i] < int(m_triangles.size()));
-
-            if (std::vector<int> neighbor_facets = neighboring_triangles(tr->children[i], second_facet_idx, second_facet_state); !neighbor_facets.empty())
-                Slic3r::append(neighbor_facets_out, std::move(neighbor_facets));
-        }
-    }
-
-    return neighbor_facets_out;
+    if (touching.second != -1)
+        process_subtriangle(touching.second);
 }
 
 void TriangleSelector::bucket_fill_select_triangles(const Vec3f& hit, int facet_start, bool propagate)
@@ -326,7 +311,23 @@ void TriangleSelector::bucket_fill_select_triangles(const Vec3f& hit, int facet_
         return;
     }
 
-    std::vector<Vec3i> all_level_neighbors = this->precompute_all_level_neighbors();
+    auto get_all_touching_triangles = [this](int facet_idx, const Vec3i &neighbors, const Vec3i &neighbors_propagated) -> std::vector<int> {
+        assert(facet_idx != -1 && facet_idx < m_triangles.size());
+        assert(this->verify_triangle_neighbors(m_triangles[facet_idx], neighbors));
+        std::vector<int> touching_triangles;
+        Vec3i            vertices = {m_triangles[facet_idx].verts_idxs[0], m_triangles[facet_idx].verts_idxs[1], m_triangles[facet_idx].verts_idxs[2]};
+        append_touching_subtriangles(neighbors(0), vertices(1), vertices(0), touching_triangles);
+        append_touching_subtriangles(neighbors(1), vertices(2), vertices(1), touching_triangles);
+        append_touching_subtriangles(neighbors(2), vertices(0), vertices(2), touching_triangles);
+
+        for (int neighbor_idx : neighbors_propagated)
+            if (neighbor_idx != -1 && !m_triangles[neighbor_idx].is_split())
+                touching_triangles.emplace_back(neighbor_idx);
+
+        return touching_triangles;
+    };
+
+    auto [neighbors, neighbors_propagated] = this->precompute_all_neighbors();
     std::vector<bool>  visited(m_triangles.size(), false);
     std::queue<int>    facet_queue;
 
@@ -338,17 +339,14 @@ void TriangleSelector::bucket_fill_select_triangles(const Vec3f& hit, int facet_
 
         if (!visited[current_facet]) {
             m_triangles[current_facet].select_by_seed_fill();
-            for (int neighbor_idx : all_level_neighbors[current_facet]) {
-                if (neighbor_idx < 0 || visited[neighbor_idx])
+
+            std::vector<int> touching_triangles = get_all_touching_triangles(current_facet, neighbors[current_facet], neighbors_propagated[current_facet]);
+            for(const int tr_idx : touching_triangles) {
+                if (tr_idx < 0 || visited[tr_idx] || m_triangles[tr_idx].get_state() != start_facet_state)
                     continue;
 
-                if (!m_triangles[neighbor_idx].is_split()) {
-                    if (m_triangles[neighbor_idx].get_state() == start_facet_state)
-                        facet_queue.push(neighbor_idx);
-                } else {
-                    for (int neighbor_facet_idx : neighboring_triangles(neighbor_idx, current_facet, start_facet_state))
-                        facet_queue.push(neighbor_facet_idx);
-                }
+                assert(!m_triangles[tr_idx].is_split());
+                facet_queue.push(tr_idx);
             }
         }
 
@@ -437,6 +435,40 @@ int TriangleSelector::neighbor_child(int itriangle, int vertexi, int vertexj, Pa
     return itriangle == -1 ? -1 : this->neighbor_child(m_triangles[itriangle], vertexi, vertexj, partition);
 }
 
+std::pair<int, int> TriangleSelector::triangle_subtriangles(int itriangle, int vertexi, int vertexj) const
+{
+    return itriangle == -1 ? std::make_pair(-1, -1) : this->triangle_subtriangles(m_triangles[itriangle], vertexi, vertexj);
+}
+
+std::pair<int, int> TriangleSelector::triangle_subtriangles(const Triangle &tr, int vertexi, int vertexj)
+{
+    if (tr.number_of_split_sides() == 0)
+        // If this triangle is not split, then there is no subtriangles touching the edge.
+        return std::make_pair(-1, -1);
+
+    // Find the triangle edge.
+    int edge = tr.verts_idxs[0] == vertexi ? 0 : tr.verts_idxs[1] == vertexi ? 1 : 2;
+    assert(tr.verts_idxs[edge] == vertexi);
+    assert(tr.verts_idxs[next_idx_modulo(edge, 3)] == vertexj);
+
+    if (tr.number_of_split_sides() == 1) {
+        return edge == next_idx_modulo(tr.special_side(), 3) ? std::make_pair(tr.children[0], tr.children[1]) :
+                                                                     std::make_pair(tr.children[edge == tr.special_side() ? 0 : 1], -1);
+    } else if (tr.number_of_split_sides() == 2) {
+        return edge == next_idx_modulo(tr.special_side(), 3) ? std::make_pair(tr.children[2], -1) :
+               edge == tr.special_side()                           ? std::make_pair(tr.children[0], tr.children[1]) :
+                                                                     std::make_pair(tr.children[2], tr.children[0]);
+    } else {
+        assert(tr.number_of_split_sides() == 3);
+        assert(tr.special_side() == 0);
+        return edge == 0 ? std::make_pair(tr.children[0], tr.children[1]) :
+               edge == 1 ? std::make_pair(tr.children[1], tr.children[2]) :
+                           std::make_pair(tr.children[2], tr.children[0]);
+    }
+
+    return std::make_pair(-1, -1);
+}
+
 // Return existing midpoint of CCW oriented side (vertexi, vertexj).
 // If itriangle == -1 or if the side sharing (vertexi, vertexj) is not split, return -1.
 int TriangleSelector::triangle_midpoint(const Triangle &tr, int vertexi, int vertexj) const
@@ -524,12 +556,8 @@ Vec3i TriangleSelector::child_neighbors(const Triangle &tr, const Vec3i &neighbo
 
     assert(child_idx >= 0 && child_idx <= tr.number_of_split_sides());
     int   i = tr.special_side();
-    int   j = i + 1;
-    if (j >= 3)
-        j = 0;
-    int   k = j + 1;
-    if (k >= 3)
-        k = 0;
+    int   j = next_idx_modulo(i, 3);
+    int   k = next_idx_modulo(j, 3);
 
     Vec3i out;
     switch (tr.number_of_split_sides()) {
@@ -612,23 +640,28 @@ Vec3i TriangleSelector::child_neighbors(const Triangle &tr, const Vec3i &neighbo
 Vec3i TriangleSelector::child_neighbors_propagated(const Triangle &tr, const Vec3i &neighbors, int child_idx) const
 {
     int i = tr.special_side();
-    int j = i + 1;
-    if (j >= 3) j = 0;
-    int k = j + 1;
-    if (k >= 3) k = 0;
+    int j = next_idx_modulo(i, 3);
+    int k = next_idx_modulo(j, 3);
 
     Vec3i out;
+    auto  replace_if_not_exists = [&out](int index_to_replace, int neighbor) {
+        if (out(index_to_replace) == -1)
+            out(index_to_replace) = neighbor;
+    };
+
     switch (tr.number_of_split_sides()) {
     case 1:
         switch (child_idx) {
         case 0:
             out(0) = neighbors(i);
-            out(1) = neighbors(j);
+            out(1) = this->neighbor_child(neighbors(j), tr.verts_idxs[k], tr.verts_idxs[j], Partition::Second);
+            replace_if_not_exists(1, neighbors(j));
             out(2) = tr.children[1];
             break;
         default:
             assert(child_idx == 1);
-            out(0) = neighbors(j);
+            out(0) = this->neighbor_child(neighbors(j), tr.verts_idxs[k], tr.verts_idxs[j], Partition::First);
+            replace_if_not_exists(0, neighbors(j));
             out(1) = neighbors(k);
             out(2) = tr.children[0];
             break;
@@ -638,20 +671,24 @@ Vec3i TriangleSelector::child_neighbors_propagated(const Triangle &tr, const Vec
     case 2:
         switch (child_idx) {
         case 0:
-            out(0) = neighbors(i);
+            out(0) = this->neighbor_child(neighbors(i), tr.verts_idxs[j], tr.verts_idxs[i], Partition::Second);
+            replace_if_not_exists(0, neighbors(i));
             out(1) = tr.children[1];
-            out(2) = neighbors(k);
+            out(2) = this->neighbor_child(neighbors(k), tr.verts_idxs[i], tr.verts_idxs[k], Partition::First);
+            replace_if_not_exists(2, neighbors(k));
             break;
         case 1:
             assert(child_idx == 1);
-            out(0) = neighbors(i);
+            out(0) = this->neighbor_child(neighbors(i), tr.verts_idxs[j], tr.verts_idxs[i], Partition::First);
+            replace_if_not_exists(0, neighbors(i));
             out(1) = tr.children[2];
             out(2) = tr.children[0];
             break;
         default:
             assert(child_idx == 2);
             out(0) = neighbors(j);
-            out(1) = neighbors(k);
+            out(1) = this->neighbor_child(neighbors(k), tr.verts_idxs[i], tr.verts_idxs[k], Partition::Second);
+            replace_if_not_exists(1, neighbors(k));
             out(2) = tr.children[1];
             break;
         }
@@ -661,18 +698,24 @@ Vec3i TriangleSelector::child_neighbors_propagated(const Triangle &tr, const Vec
         assert(tr.special_side() == 0);
         switch (child_idx) {
         case 0:
-            out(0) = neighbors(0);
+            out(0) = this->neighbor_child(neighbors(0), tr.verts_idxs[1], tr.verts_idxs[0], Partition::Second);
+            replace_if_not_exists(0, neighbors(0));
             out(1) = tr.children[3];
-            out(2) = neighbors(2);
+            out(2) = this->neighbor_child(neighbors(2), tr.verts_idxs[0], tr.verts_idxs[2], Partition::First);
+            replace_if_not_exists(2, neighbors(2));
             break;
         case 1:
-            out(0) = neighbors(0);
-            out(1) = neighbors(1);
+            out(0) = this->neighbor_child(neighbors(0), tr.verts_idxs[1], tr.verts_idxs[0], Partition::First);
+            replace_if_not_exists(0, neighbors(0));
+            out(1) = this->neighbor_child(neighbors(1), tr.verts_idxs[2], tr.verts_idxs[1], Partition::Second);
+            replace_if_not_exists(1, neighbors(1));
             out(2) = tr.children[3];
             break;
         case 2:
-            out(0) = neighbors(1);
-            out(1) = neighbors(2);
+            out(0) = this->neighbor_child(neighbors(1), tr.verts_idxs[2], tr.verts_idxs[1], Partition::First);
+            replace_if_not_exists(0, neighbors(1));
+            out(1) = this->neighbor_child(neighbors(2), tr.verts_idxs[0], tr.verts_idxs[2], Partition::Second);
+            replace_if_not_exists(1, neighbors(2));
             out(2) = tr.children[3];
             break;
         default:
@@ -886,13 +929,13 @@ void TriangleSelector::undivide_triangle(int facet_idx)
     Triangle& tr = m_triangles[facet_idx];
 
     if (tr.is_split()) {
-        for (int i=0; i<=tr.number_of_split_sides(); ++i) {
+        for (int i = 0; i <= tr.number_of_split_sides(); ++i) {
             int       child    = tr.children[i];
             Triangle &child_tr = m_triangles[child];
             assert(child_tr.valid());
             undivide_triangle(child);
-            for (int i = 0; i < 3; ++ i) {
-                int     iv = child_tr.verts_idxs[i];
+            for (int j = 0; j < 3; ++j) {
+                int     iv = child_tr.verts_idxs[j];
                 Vertex &v  = m_vertices[iv];
                 assert(v.ref_cnt > 0);
                 if (-- v.ref_cnt == 0) {
@@ -1231,7 +1274,7 @@ void TriangleSelector::get_facets_strict_recursive(
         this->get_facets_split_by_tjoints({tr.verts_idxs[0], tr.verts_idxs[1], tr.verts_idxs[2]}, neighbors, out_triangles);
 }
 
-void TriangleSelector::get_facets_split_by_tjoints(const Vec3i vertices, const Vec3i neighbors, std::vector<stl_triangle_vertex_indices> &out_triangles) const
+void TriangleSelector::get_facets_split_by_tjoints(const Vec3i &vertices, const Vec3i &neighbors, std::vector<stl_triangle_vertex_indices> &out_triangles) const
 {
 // Export this triangle, but first collect the T-joint vertices along its edges.
     Vec3i midpoints(
@@ -1393,9 +1436,10 @@ std::pair<std::vector<std::pair<int, int>>, std::vector<bool>> TriangleSelector:
     return out.data;
 }
 
-void TriangleSelector::deserialize(const std::pair<std::vector<std::pair<int, int>>, std::vector<bool>> &data)
+void TriangleSelector::deserialize(const std::pair<std::vector<std::pair<int, int>>, std::vector<bool>> &data, bool needs_reset)
 {
-    reset(); // dump any current state
+    if (needs_reset)
+        reset(); // dump any current state
 
     // Reserve number of triangles as if each triangle was saved with 4 bits.
     // With MMU painting this estimate may be somehow low, but better than nothing.
diff --git a/src/libslic3r/TriangleSelector.hpp b/src/libslic3r/TriangleSelector.hpp
index 643daba45..eeb479dee 100644
--- a/src/libslic3r/TriangleSelector.hpp
+++ b/src/libslic3r/TriangleSelector.hpp
@@ -22,8 +22,8 @@ public:
         POINTER
     };
 
-    [[nodiscard]] std::vector<Vec3i> precompute_all_level_neighbors() const;
-    void precompute_all_level_neighbors_recursive(const int facet_idx, const Vec3i &neighbors, const Vec3i &neighbors_propagated, std::vector<Vec3i> &neighbors_out) const;
+    std::pair<std::vector<Vec3i>, std::vector<Vec3i>> precompute_all_neighbors() const;
+    void precompute_all_neighbors_recursive(int facet_idx, const Vec3i &neighbors, const Vec3i &neighbors_propagated, std::vector<Vec3i> &neighbors_out, std::vector<Vec3i> &neighbors_normal_out) const;
 
     // Set a limit to the edge length, below which the edge will not be split by select_patch().
     // Called by select_patch() internally. Made public for debugging purposes, see TriangleSelectorGUI::render_debug().
@@ -37,10 +37,6 @@ public:
     [[nodiscard]] int select_unsplit_triangle(const Vec3f &hit, int facet_idx) const;
     [[nodiscard]] int select_unsplit_triangle(const Vec3f &hit, int facet_idx, const Vec3i &neighbors) const;
 
-    [[nodiscard]] bool are_triangles_touching(int first_facet_idx, int second_facet_idx) const;
-
-    [[nodiscard]] std::vector<int> neighboring_triangles(int first_facet_idx, int second_facet_idx, EnforcerBlockerType second_facet_state) const;
-
     // Select all triangles fully inside the circle, subdivide where needed.
     void select_patch(const Vec3f        &hit,         // point where to start
                       int                 facet_start, // facet of the original mesh (unsplit) that the hit point belongs to
@@ -60,7 +56,7 @@ public:
                                     bool         propagate);        // if bucket fill is propagated to neighbor faces or if it fills the only facet of the modified mesh that the hit point belongs to.
 
     bool                 has_facets(EnforcerBlockerType state) const;
-    static bool          has_facets(const std::pair<std::vector<std::pair<int, int>>, std::vector<bool>> &data, const EnforcerBlockerType test_state);
+    static bool          has_facets(const std::pair<std::vector<std::pair<int, int>>, std::vector<bool>> &data, EnforcerBlockerType test_state);
     int                  num_facets(EnforcerBlockerType state) const;
     // Get facets at a given state. Don't triangulate T-joints.
     indexed_triangle_set get_facets(EnforcerBlockerType state) const;
@@ -81,7 +77,7 @@ public:
     std::pair<std::vector<std::pair<int, int>>, std::vector<bool>> serialize() const;
 
     // Load serialized data. Assumes that correct mesh is loaded.
-    void deserialize(const std::pair<std::vector<std::pair<int, int>>, std::vector<bool>> &data);
+    void deserialize(const std::pair<std::vector<std::pair<int, int>>, std::vector<bool>> &data, bool needs_reset = true);
 
     // For all triangles, remove the flag indicating that the triangle was selected by seed fill.
     void seed_fill_unselect_all_triangles();
@@ -128,11 +124,11 @@ protected:
         bool is_selected_by_seed_fill() const { assert(! is_split()); return m_selected_by_seed_fill; }
 
         // Is this triangle valid or marked to be removed?
-        bool valid() const throw() { return m_valid; }
+        bool valid() const noexcept { return m_valid; }
         // Get info on how it's split.
-        bool is_split() const throw() { return number_of_split_sides() != 0; }
-        int number_of_split_sides() const throw() { return number_of_splits; }
-        int special_side() const throw() { assert(is_split()); return special_side_idx; }
+        bool is_split() const noexcept { return number_of_split_sides() != 0; }
+        int number_of_split_sides() const noexcept { return number_of_splits; }
+        int special_side() const noexcept { assert(is_split()); return special_side_idx; }
 
     private:
         friend TriangleSelector;
@@ -205,7 +201,7 @@ private:
     void remove_useless_children(int facet_idx); // No hidden meaning. Triangles are meant.
     bool is_pointer_in_triangle(int facet_idx) const;
     bool is_edge_inside_cursor(int facet_idx) const;
-    int  push_triangle(int a, int b, int c, int source_triangle, const EnforcerBlockerType state = EnforcerBlockerType{0});
+    int  push_triangle(int a, int b, int c, int source_triangle, EnforcerBlockerType state = EnforcerBlockerType{0});
     void perform_split(int facet_idx, const Vec3i &neighbors, EnforcerBlockerType old_state);
     Vec3i child_neighbors(const Triangle &tr, const Vec3i &neighbors, int child_idx) const;
     Vec3i child_neighbors_propagated(const Triangle &tr, const Vec3i &neighbors, int child_idx) const;
@@ -221,6 +217,11 @@ private:
     int triangle_midpoint(int itriangle, int vertexi, int vertexj) const;
     int triangle_midpoint_or_allocate(int itriangle, int vertexi, int vertexj);
 
+    static std::pair<int, int> triangle_subtriangles(const Triangle &tr, int vertexi, int vertexj);
+    std::pair<int, int>        triangle_subtriangles(int itriangle, int vertexi, int vertexj) const;
+
+    void append_touching_subtriangles(int itriangle, int vertexi, int vertexj, std::vector<int> &touching_subtriangles_out) const;
+
 #ifndef NDEBUG
     bool verify_triangle_neighbors(const Triangle& tr, const Vec3i& neighbors) const;
     bool verify_triangle_midpoints(const Triangle& tr) const;
@@ -231,7 +232,7 @@ private:
         const Vec3i                                 &neighbors,
         EnforcerBlockerType                          state,
         std::vector<stl_triangle_vertex_indices>    &out_triangles) const;
-    void get_facets_split_by_tjoints(const Vec3i vertices, const Vec3i neighbors, std::vector<stl_triangle_vertex_indices> &out_triangles) const;
+    void get_facets_split_by_tjoints(const Vec3i &vertices, const Vec3i &neighbors, std::vector<stl_triangle_vertex_indices> &out_triangles) const;
 
     int m_free_triangles_head { -1 };
     int m_free_vertices_head { -1 };
diff --git a/src/slic3r/GUI/GLCanvas3D.cpp b/src/slic3r/GUI/GLCanvas3D.cpp
index 6a40a4e44..6f34f4051 100644
--- a/src/slic3r/GUI/GLCanvas3D.cpp
+++ b/src/slic3r/GUI/GLCanvas3D.cpp
@@ -1407,10 +1407,6 @@ void GLCanvas3D::render()
     if (!is_initialized() && !init())
         return;
 
-#if ENABLE_RENDER_STATISTICS
-    auto start_time = std::chrono::high_resolution_clock::now();
-#endif // ENABLE_RENDER_STATISTICS
-
     if (wxGetApp().plater()->get_bed().get_shape().empty()) {
         // this happens at startup when no data is still saved under <>\AppData\Roaming\Slic3rPE
         post_event(SimpleEvent(EVT_GLCANVAS_UPDATE_BED_SHAPE));
@@ -1505,19 +1501,12 @@ void GLCanvas3D::render()
     // draw overlays
     _render_overlays();
 
-#if ENABLE_RENDER_STATISTICS
     if (wxGetApp().plater()->is_render_statistic_dialog_visible()) {
         ImGuiWrapper& imgui = *wxGetApp().imgui();
         imgui.begin(std::string("Render statistics"), ImGuiWindowFlags_AlwaysAutoResize | ImGuiWindowFlags_NoResize | ImGuiWindowFlags_NoCollapse);
-        imgui.text("Last frame:");
+        imgui.text("FPS (SwapBuffers() calls per second):");
         ImGui::SameLine();
-        int64_t average = m_render_stats.get_average();
-        imgui.text(std::to_string(average));
-        ImGui::SameLine();
-        imgui.text("ms");
-        imgui.text("FPS:");
-        ImGui::SameLine();
-        imgui.text(std::to_string((average == 0) ? 0 : static_cast<int>(1000.0f / static_cast<float>(average))));
+        imgui.text(std::to_string(m_render_stats.get_fps_and_reset_if_needed()));
         ImGui::Separator();
         imgui.text("Compressed textures:");
         ImGui::SameLine();
@@ -1527,7 +1516,6 @@ void GLCanvas3D::render()
         imgui.text(std::to_string(OpenGLManager::get_gl_info().get_max_tex_size()));
         imgui.end();
     }
-#endif // ENABLE_RENDER_STATISTICS
 
 #if ENABLE_PROJECT_DIRTY_STATE_DEBUG_WINDOW
     if (wxGetApp().is_editor() && wxGetApp().plater()->is_view3D_shown())
@@ -1574,11 +1562,7 @@ void GLCanvas3D::render()
     wxGetApp().imgui()->render();
 
     m_canvas->SwapBuffers();
-
-#if ENABLE_RENDER_STATISTICS
-    auto end_time = std::chrono::high_resolution_clock::now();
-    m_render_stats.add_frame(std::chrono::duration_cast<std::chrono::milliseconds>(end_time - start_time).count());
-#endif // ENABLE_RENDER_STATISTICS
+    m_render_stats.increment_fps_counter();
 }
 
 void GLCanvas3D::render_thumbnail(ThumbnailData& thumbnail_data, unsigned int w, unsigned int h, const ThumbnailsParams& thumbnail_params, Camera::EType camera_type)
@@ -2592,15 +2576,11 @@ void GLCanvas3D::on_key(wxKeyEvent& evt)
     {
         if (!m_gizmos.on_key(evt)) {
             if (evt.GetEventType() == wxEVT_KEY_UP) {
-#if ENABLE_RENDER_STATISTICS
                 if (evt.ShiftDown() && evt.ControlDown() && keyCode == WXK_SPACE) {
                     wxGetApp().plater()->toggle_render_statistic_dialog();
                     m_dirty = true;
                 }
                 if (m_tab_down && keyCode == WXK_TAB && !evt.HasAnyModifiers()) {
-#else
-                if (m_tab_down && keyCode == WXK_TAB && !evt.HasAnyModifiers()) {
-#endif // ENABLE_RENDER_STATISTICS
                     // Enable switching between 3D and Preview with Tab
                     // m_canvas->HandleAsNavigationKey(evt);   // XXX: Doesn't work in some cases / on Linux
                     post_event(SimpleEvent(EVT_GLCANVAS_TAB));
@@ -3437,6 +3417,7 @@ void GLCanvas3D::do_move(const std::string& snapshot_type)
             m_selection.translate(i.first, i.second, shift);
             m->translate_instance(i.second, shift);
         }
+        wxGetApp().obj_list()->update_info_items(static_cast<size_t>(i.first));
     }
 
     // if the selection is not valid to allow for layer editing after the move, we need to turn off the tool if it is running
@@ -3517,6 +3498,7 @@ void GLCanvas3D::do_rotate(const std::string& snapshot_type)
             m_selection.translate(i.first, i.second, shift);
             m->translate_instance(i.second, shift);
         }
+        wxGetApp().obj_list()->update_info_items(static_cast<size_t>(i.first));
     }
 
     if (!done.empty())
@@ -3584,6 +3566,7 @@ void GLCanvas3D::do_scale(const std::string& snapshot_type)
             m_selection.translate(i.first, i.second, shift);
             m->translate_instance(i.second, shift);
         }
+        wxGetApp().obj_list()->update_info_items(static_cast<size_t>(i.first));
     }
 
     if (!done.empty())
diff --git a/src/slic3r/GUI/GLCanvas3D.hpp b/src/slic3r/GUI/GLCanvas3D.hpp
index b46778c39..b5fa86235 100644
--- a/src/slic3r/GUI/GLCanvas3D.hpp
+++ b/src/slic3r/GUI/GLCanvas3D.hpp
@@ -305,25 +305,27 @@ class GLCanvas3D
         ObjectClashed
     };
 
-#if ENABLE_RENDER_STATISTICS
     class RenderStats
     {
-        std::queue<std::pair<int64_t, int64_t>> m_frames;
-        int64_t m_curr_total{ 0 };
-
+    private:
+        std::chrono::time_point<std::chrono::high_resolution_clock> m_measuring_start;
+        int m_fps_out = -1;
+        int m_fps_running = 0;
     public:
-        void add_frame(int64_t frame) {
-            int64_t now = GLCanvas3D::timestamp_now();
-            if (!m_frames.empty() && now - m_frames.front().first > 1000) {
-                m_curr_total -= m_frames.front().second;
-                m_frames.pop();
+        void increment_fps_counter() { ++m_fps_running; }
+        int get_fps() { return m_fps_out; }
+        int get_fps_and_reset_if_needed() {
+            auto cur_time = std::chrono::high_resolution_clock::now();
+            int elapsed_ms = std::chrono::duration_cast<std::chrono::milliseconds>(cur_time-m_measuring_start).count();
+            if (elapsed_ms > 1000  || m_fps_out == -1) {
+                m_measuring_start = cur_time;
+                m_fps_out = int (1000. * m_fps_running / elapsed_ms);
+                m_fps_running = 0;
             }
-            m_curr_total += frame;
-            m_frames.push({ now, frame });
+            return m_fps_out;
         }
-        int64_t get_average() const { return m_frames.empty() ? 0 : m_curr_total / m_frames.size(); }
+
     };
-#endif // ENABLE_RENDER_STATISTICS
 
     class Labels
     {
@@ -455,9 +457,7 @@ private:
     bool m_show_picking_texture;
 #endif // ENABLE_RENDER_PICKING_PASS
 
-#if ENABLE_RENDER_STATISTICS
     RenderStats m_render_stats;
-#endif // ENABLE_RENDER_STATISTICS
 
     int m_imgui_undo_redo_hovered_pos{ -1 };
     int m_mouse_wheel{ 0 };
diff --git a/src/slic3r/GUI/GUI_ObjectList.cpp b/src/slic3r/GUI/GUI_ObjectList.cpp
index a013b8dd2..834e45898 100644
--- a/src/slic3r/GUI/GUI_ObjectList.cpp
+++ b/src/slic3r/GUI/GUI_ObjectList.cpp
@@ -2386,15 +2386,28 @@ void ObjectList::part_selection_changed()
 
                 if (type == itInfo) {
                     InfoItemType info_type = m_objects_model->GetInfoItemType(item);
-                    if (info_type != InfoItemType::VariableLayerHeight) {
+                    switch (info_type)
+                    {
+                    case InfoItemType::VariableLayerHeight:
+                    {
+                        wxGetApp().plater()->toggle_layers_editing(true);
+                        break;
+                    }
+                    case InfoItemType::CustomSupports:
+                    case InfoItemType::CustomSeam:
+                    case InfoItemType::MmuSegmentation:
+                    {
                         GLGizmosManager::EType gizmo_type = info_type == InfoItemType::CustomSupports ? GLGizmosManager::EType::FdmSupports :
-                                                            info_type == InfoItemType::CustomSeam     ? GLGizmosManager::EType::Seam :
-                                                                                                        GLGizmosManager::EType::MmuSegmentation;
+                                                            info_type == InfoItemType::CustomSeam ? GLGizmosManager::EType::Seam :
+                                                            GLGizmosManager::EType::MmuSegmentation;
                         GLGizmosManager& gizmos_mgr = wxGetApp().plater()->canvas3D()->get_gizmos_manager();
                         if (gizmos_mgr.get_current_type() != gizmo_type)
                             gizmos_mgr.open_gizmo(gizmo_type);
-                    } else
-                        wxGetApp().plater()->toggle_layers_editing(true);
+                        break;
+                    }
+                    case InfoItemType::Sinking: { break; }
+                    default: { break; }
+                    }
                 }
             }
             else {
@@ -2520,6 +2533,7 @@ void ObjectList::update_info_items(size_t obj_idx)
     for (InfoItemType type : {InfoItemType::CustomSupports,
                               InfoItemType::CustomSeam,
                               InfoItemType::MmuSegmentation,
+                              InfoItemType::Sinking,
                               InfoItemType::VariableLayerHeight}) {
         wxDataViewItem item = m_objects_model->GetInfoItemByType(item_obj, type);
         bool shows = item.IsOk();
@@ -2542,6 +2556,13 @@ void ObjectList::update_info_items(size_t obj_idx)
             should_show = printer_technology() == ptFFF
                        && ! model_object->layer_height_profile.empty();
             break;
+        case InfoItemType::Sinking:
+        {
+            const BoundingBoxf3& box = model_object->bounding_box();
+            should_show = printer_technology() == ptFFF &&
+                        box.min.z() < SINKING_Z_THRESHOLD && box.max.z() > SINKING_Z_THRESHOLD;
+            break;
+        }
         default: break;
         }
 
diff --git a/src/slic3r/GUI/Gizmos/GLGizmoFdmSupports.cpp b/src/slic3r/GUI/Gizmos/GLGizmoFdmSupports.cpp
index f540e30aa..68f0f3f99 100644
--- a/src/slic3r/GUI/Gizmos/GLGizmoFdmSupports.cpp
+++ b/src/slic3r/GUI/Gizmos/GLGizmoFdmSupports.cpp
@@ -346,7 +346,8 @@ void GLGizmoFdmSupports::update_from_model_object()
         const TriangleMesh* mesh = &mv->mesh();
 
         m_triangle_selectors.emplace_back(std::make_unique<TriangleSelectorGUI>(*mesh));
-        m_triangle_selectors.back()->deserialize(mv->supported_facets.get_data());
+        // Reset of TriangleSelector is done inside TriangleSelectorGUI's constructor, so we don't need it to perform it again in deserialize().
+        m_triangle_selectors.back()->deserialize(mv->supported_facets.get_data(), false);
         m_triangle_selectors.back()->request_update_render_data();
     }
 }
diff --git a/src/slic3r/GUI/Gizmos/GLGizmoMmuSegmentation.cpp b/src/slic3r/GUI/Gizmos/GLGizmoMmuSegmentation.cpp
index a35d8c071..f6e7708fa 100644
--- a/src/slic3r/GUI/Gizmos/GLGizmoMmuSegmentation.cpp
+++ b/src/slic3r/GUI/Gizmos/GLGizmoMmuSegmentation.cpp
@@ -544,7 +544,8 @@ void GLGizmoMmuSegmentation::init_model_triangle_selectors()
 
         int extruder_idx = (mv->extruder_id() > 0) ? mv->extruder_id() - 1 : 0;
         m_triangle_selectors.emplace_back(std::make_unique<TriangleSelectorMmGui>(*mesh, m_modified_extruders_colors, m_original_extruders_colors[size_t(extruder_idx)]));
-        m_triangle_selectors.back()->deserialize(mv->mmu_segmentation_facets.get_data());
+        // Reset of TriangleSelector is done inside TriangleSelectorMmGUI's constructor, so we don't need it to perform it again in deserialize().
+        m_triangle_selectors.back()->deserialize(mv->mmu_segmentation_facets.get_data(), false);
         m_triangle_selectors.back()->request_update_render_data();
     }
     m_original_volumes_extruder_idxs = get_extruder_id_for_volumes(*mo);
diff --git a/src/slic3r/GUI/Gizmos/GLGizmoSeam.cpp b/src/slic3r/GUI/Gizmos/GLGizmoSeam.cpp
index 5f0c6b52d..8854c1a7e 100644
--- a/src/slic3r/GUI/Gizmos/GLGizmoSeam.cpp
+++ b/src/slic3r/GUI/Gizmos/GLGizmoSeam.cpp
@@ -256,7 +256,8 @@ void GLGizmoSeam::update_from_model_object()
         const TriangleMesh* mesh = &mv->mesh();
 
         m_triangle_selectors.emplace_back(std::make_unique<TriangleSelectorGUI>(*mesh));
-        m_triangle_selectors.back()->deserialize(mv->seam_facets.get_data());
+        // Reset of TriangleSelector is done inside TriangleSelectorGUI's constructor, so we don't need it to perform it again in deserialize().
+        m_triangle_selectors.back()->deserialize(mv->seam_facets.get_data(), false);
         m_triangle_selectors.back()->request_update_render_data();
     }
 }
diff --git a/src/slic3r/GUI/ObjectDataViewModel.cpp b/src/slic3r/GUI/ObjectDataViewModel.cpp
index 3eb0cd5c9..9f48bcc3c 100644
--- a/src/slic3r/GUI/ObjectDataViewModel.cpp
+++ b/src/slic3r/GUI/ObjectDataViewModel.cpp
@@ -65,6 +65,7 @@ ObjectDataViewModelNode::ObjectDataViewModelNode(ObjectDataViewModelNode* parent
     m_name           = info_type == InfoItemType::CustomSupports  ? _L("Paint-on supports") :
                        info_type == InfoItemType::CustomSeam      ? _L("Paint-on seam") :
                        info_type == InfoItemType::MmuSegmentation ? _L("Paint-on segmentation") :
+                       info_type == InfoItemType::Sinking         ? _L("Sinking") :
                                                                     _L("Variable layer height");
     m_info_item_type = info_type;
 }
diff --git a/src/slic3r/GUI/ObjectDataViewModel.hpp b/src/slic3r/GUI/ObjectDataViewModel.hpp
index da251ef84..1cf10faf5 100644
--- a/src/slic3r/GUI/ObjectDataViewModel.hpp
+++ b/src/slic3r/GUI/ObjectDataViewModel.hpp
@@ -51,6 +51,7 @@ enum class InfoItemType
     CustomSupports,
     CustomSeam,
     MmuSegmentation,
+    Sinking,
     VariableLayerHeight
 };
 
diff --git a/src/slic3r/GUI/Plater.cpp b/src/slic3r/GUI/Plater.cpp
index 649ad88cf..9b59026dc 100644
--- a/src/slic3r/GUI/Plater.cpp
+++ b/src/slic3r/GUI/Plater.cpp
@@ -1565,9 +1565,7 @@ struct Plater::priv
     std::string                 label_btn_export;
     std::string                 label_btn_send;
 
-#if ENABLE_RENDER_STATISTICS
     bool                        show_render_statistic_dialog{ false };
-#endif // ENABLE_RENDER_STATISTICS
 
     static const std::regex pattern_bundle;
     static const std::regex pattern_3mf;
@@ -6516,7 +6514,6 @@ void Plater::enter_gizmos_stack() { p->enter_gizmos_stack(); }
 void Plater::leave_gizmos_stack() { p->leave_gizmos_stack(); }
 bool Plater::inside_snapshot_capture() { return p->inside_snapshot_capture(); }
 
-#if ENABLE_RENDER_STATISTICS
 void Plater::toggle_render_statistic_dialog()
 {
     p->show_render_statistic_dialog = !p->show_render_statistic_dialog;
@@ -6526,7 +6523,6 @@ bool Plater::is_render_statistic_dialog_visible() const
 {
     return p->show_render_statistic_dialog;
 }
-#endif // ENABLE_RENDER_STATISTICS
 
 // Wrapper around wxWindow::PopupMenu to suppress error messages popping out while tracking the popup menu.
 bool Plater::PopupMenu(wxMenu *menu, const wxPoint& pos)
diff --git a/src/slic3r/GUI/Plater.hpp b/src/slic3r/GUI/Plater.hpp
index 0b499725e..9c37dfe2b 100644
--- a/src/slic3r/GUI/Plater.hpp
+++ b/src/slic3r/GUI/Plater.hpp
@@ -404,10 +404,8 @@ public:
 
     bool inside_snapshot_capture();
 
-#if ENABLE_RENDER_STATISTICS
     void toggle_render_statistic_dialog();
     bool is_render_statistic_dialog_visible() const;
-#endif // ENABLE_RENDER_STATISTICS
 
 	// Wrapper around wxWindow::PopupMenu to suppress error messages popping out while tracking the popup menu.
 	bool PopupMenu(wxMenu *menu, const wxPoint& pos = wxDefaultPosition);
diff --git a/src/slic3r/GUI/Selection.cpp b/src/slic3r/GUI/Selection.cpp
index e8d7cc621..ce44e5442 100644
--- a/src/slic3r/GUI/Selection.cpp
+++ b/src/slic3r/GUI/Selection.cpp
@@ -841,8 +841,12 @@ void Selection::scale(const Vec3d& scale, TransformationType transformation_type
 
     for (unsigned int i : m_list) {
         GLVolume &v = *(*m_volumes)[i];
-        if (!is_sla)
-            is_any_volume_sinking |= !v.is_modifier && std::find(m_cache.sinking_volumes.begin(), m_cache.sinking_volumes.end(), i) != m_cache.sinking_volumes.end();
+        if (!is_sla) {
+            if (v.is_modifier)
+                is_any_volume_sinking = true;
+            else
+                is_any_volume_sinking |= std::find(m_cache.sinking_volumes.begin(), m_cache.sinking_volumes.end(), i) != m_cache.sinking_volumes.end();
+        }
         if (is_single_full_instance()) {
             if (transformation_type.relative()) {
                 Transform3d m = Geometry::assemble_transform(Vec3d::Zero(), Vec3d::Zero(), scale);
@@ -2123,7 +2127,7 @@ void Selection::ensure_on_bed()
 
     for (GLVolume* volume : *m_volumes) {
         if (!volume->is_wipe_tower && !volume->is_modifier) {
-            double min_z = volume->transformed_convex_hull_bounding_box().min(2);
+            const double min_z = volume->transformed_convex_hull_bounding_box().min.z();
             std::pair<int, int> instance = std::make_pair(volume->object_idx(), volume->instance_idx());
             InstancesToZMap::iterator it = instances_min_z.find(instance);
             if (it == instances_min_z.end())
diff --git a/src/slic3r/GUI/Tab.cpp b/src/slic3r/GUI/Tab.cpp
index a91b17dd3..989ff04c4 100644
--- a/src/slic3r/GUI/Tab.cpp
+++ b/src/slic3r/GUI/Tab.cpp
@@ -3059,6 +3059,7 @@ void Tab::load_current_preset()
                         if (!wxGetApp().tabs_as_menu()) {
                             std::string bmp_name = tab->type() == Slic3r::Preset::TYPE_FILAMENT      ? "spool" :
                                                    tab->type() == Slic3r::Preset::TYPE_SLA_MATERIAL  ? "resin" : "cog";
+                            tab->Hide(); // #ys_WORKAROUND : Hide tab before inserting to avoid unwanted rendering of the tab
                             dynamic_cast<Notebook*>(wxGetApp().tab_panel())->InsertPage(wxGetApp().tab_panel()->FindPage(this), tab, tab->title(), bmp_name);
                         }
                         else