1....Binary Threshold...

% Main script to read an image, convert it to binary, and display the result
% altime using gray image
img = imread('ccc.jpg');
% Call convert2binary() function to convert the image to binary
binary_image = convert2binary(img);
% Function to convert an image to binary based on thresholding
function [binary] = convert2binary(img)
    [x, y, z] = size(img);
    % If the input is an RGB image, convert it to grayscale
    if z == 3
        img = 0.2989 * img(:,:,1) + 0.5870 * img(:,:,2) + 0.1140 * img(:,:,3); % Manual grayscale conversion
    end
    % Convert image class from 'uint8' to 'double'
    img = double(img);
    % Calculate the sum of all gray level pixel values
    total_sum = sum(img(:));
    % Calculate threshold as the average pixel value
    threshold = total_sum / (x * y);
    % Initialize the binary image matrix
    binary = zeros(x, y);
    % Apply thresholding
    for i = 1:x
        for j = 1:y
            if img(i, j) >= threshold
                binary(i, j) = 1;
            else
                binary(i, j) = 0;
            end
        end
    end
end
figure;
subplot(1, 2, 1);
imshow(img);
title('Original Image');
subplot(1, 2, 2);
imshow(binary_image);
title('Binary Image');

2....Contrast...

% Load the image
img = imread('color.jpg');  % Replace with the name of your image file

% Convert the image to grayscale (if not already)
gray_img = rgb2gray(img);

% Get the minimum and maximum pixel values of the grayscale image
r_min = double(min(gray_img(:)));  % Minimum pixel value
r_max = double(max(gray_img(:)));  % Maximum pixel value

% Set the desired output range
s_min = 0;  % Minimum output pixel value
s_max = 255;  % Maximum output pixel value

% Apply contrast stretching
contrast_stretched_img = uint8((double(gray_img) - r_min) / (r_max - r_min) * (s_max - s_min) + s_min);

% Display the original and contrast stretched images
figure;
subplot(1, 2, 1);
imshow(gray_img);
title('Original Image');

subplot(1, 2, 2);
imshow(contrast_stretched_img);
title('Contrast Stretched Image');

3.....Feature extraction...

clc; clear; close all;

% Step 1: Read the color image
image = imread('apple.jpeg');  % Replace with your image file path

% Step 2: Extract RGB components
red_component = image(:,:,1);    % Red component (1st channel)
green_component = image(:,:,2);  % Green component (2nd channel)
blue_component = image(:,:,3);   % Blue component (3rd channel)

% Step 3: Manually convert to grayscale using the formula
gray_image_manual = 0.2989 * double(red_component) + 0.5870 * double(green_component) + 0.1140 * double(blue_component);
gray_image_manual = uint8(gray_image_manual);  % Convert back to uint8 for display

% Step 4: Display the original and component images
figure;

% Original Image
subplot(2,3,1);
imshow(image);
title('Original Color Image');

% Red Component
subplot(2,3,2);
imshow(red_component);
title('Red Component');

% Green Component
subplot(2,3,3);
imshow(green_component);
title('Green Component');

% Blue Component
subplot(2,3,4);
imshow(blue_component);
title('Blue Component');

% Grayscale Image (Manually Converted)
subplot(2,3,5);
imshow(gray_image_manual);
title('Grayscale Image (Manual Conversion)');

4....Histogram....

% Load the image
image = imread('lab practise/apple.jpeg'); % Replace 'your_image.jpg' with the actual image path

% Check if the image is grayscale or RGB
if size(image, 3) == 1
    % For Grayscale image
    figure;
    histogram(image, 256); % 256 bins for pixel intensity
    title('Histogram of Grayscale Image');
    xlabel('Pixel Intensity');
    ylabel('Frequency');
else
    % For RGB image
    redChannel = image(:,:,1); % Red channel
    greenChannel = image(:,:,2); % Green channel
    blueChannel = image(:,:,3); % Blue channel

    % Plot histograms for each channel
    figure;
    subplot(3, 1, 1);
    histogram(redChannel, 256, 'FaceColor', 'r');
    title('Red Channel Histogram');
    xlabel('Pixel Intensity');
    ylabel('Frequency');

    subplot(3, 1, 2);
    histogram(greenChannel, 256, 'FaceColor', 'g');
    title('Green Channel Histogram');
    xlabel('Pixel Intensity');
    ylabel('Frequency');

    subplot(3, 1, 3);
    histogram(blueChannel, 256, 'FaceColor', 'b');
    title('Blue Channel Histogram');
    xlabel('Pixel Intensity');
    ylabel('Frequency');
end


5....Negative...

% Load the image
img = imread('flower.jpg'); % Replace 'negativeim.jpeg' with your image file name

% Convert the image to double precision to prevent data loss during inversion
img = double(img);

% Maximum pixel value for an 8-bit image
L = 255;

% Generate the negative image by applying the equation s = L - 1 - r
negative_img = L - 1 - img;

% Convert back to uint8 for display
negative_img = uint8(negative_img);

% Display the original and negative images
figure;
subplot(1, 2, 1);
imshow(uint8(img)); % Display the original image
title('Original Image');

subplot(1, 2, 2);
imshow(negative_img); % Display the negative image
title('Negative Image');

6....rgb2gray...

% Read the RGB image
image = imread('color.jpg');

% Separate the color channels
R = image(:, :, 1);
G = image(:, :, 2);
B = image(:, :, 3);

% Convert RGB to grayscale
grayImage = uint8(0.2989 * double(R) + 0.5870 * double(G) + 0.1140 * double(B));

% Convert grayscale back to RGB by replicating grayscale values across R, G, and B channels
rgbImage = cat(3, R, G, B);

% Display images
figure;
subplot(2, 3, 1), imshow(image), title('Original RGB Image');
subplot(2, 3, 2), imshow(R), title('Red Channel');
subplot(2, 3, 3), imshow(G), title('Green Channel');
subplot(2, 3, 4), imshow(B), title('Blue Channel');
subplot(2, 3, 5), imshow(grayImage), title('Grayscale Image');
subplot(2, 3, 6), imshow(rgbImage), title('Reconstructed RGB Image');

....Robert+sobel+prewitt....
% Read and convert the input image to grayscale
inputImage = imread('ccc.jpg'); % Replace with your image file
grayImage = rgb2gray(inputImage);       % Convert to grayscale if it's a color image
grayImage = double(grayImage);          % Convert to double for computation

% Initialize the size of the image
[rows, cols] = size(grayImage);

% Define the kernels
RobertX = [1 0; 0 -1];  % Robert Cross gradient kernel (X direction)
RobertY = [0 1; -1 0];  % Robert Cross gradient kernel (Y direction)

PrewittX = [-1 0 1; -1 0 1; -1 0 1];  % Prewitt kernel (X direction)
PrewittY = [-1 -1 -1; 0 0 0; 1 1 1];  % Prewitt kernel (Y direction)

SobelX = [-1 0 1; -2 0 2; -1 0 1];    % Sobel kernel (X direction)
SobelY = [-1 -2 -1; 0 0 0; 1 2 1];    % Sobel kernel (Y direction)

% Initialize the gradient magnitude images
RobertEdge = zeros(rows, cols);
PrewittEdge = zeros(rows, cols);
SobelEdge = zeros(rows, cols);

% Apply the Robert Operator (X and Y gradients)
for i = 1:rows-1
    for j = 1:cols-1
        % Apply Robert X kernel
        GxR = sum(sum(grayImage(i:i+1, j:j+1) .* RobertX));
        % Apply Robert Y kernel
        GyR = sum(sum(grayImage(i:i+1, j:j+1) .* RobertY));
        % Compute gradient magnitude
        RobertEdge(i,j) = sqrt(GxR^2 + GyR^2);
    end
end

% Apply the Prewitt Operator (X and Y gradients)
for i = 2:rows-1
    for j = 2:cols-1
        % Apply Prewitt X kernel
        GxP = sum(sum(grayImage(i-1:i+1, j-1:j+1) .* PrewittX));
        % Apply Prewitt Y kernel
        GyP = sum(sum(grayImage(i-1:i+1, j-1:j+1) .* PrewittY));
        % Compute gradient magnitude
        PrewittEdge(i,j) = sqrt(GxP^2 + GyP^2);
    end
end

% Apply the Sobel Operator (X and Y gradients)
for i = 2:rows-1
    for j = 2:cols-1
        % Apply Sobel X kernel
        GxS = sum(sum(grayImage(i-1:i+1, j-1:j+1) .* SobelX));
        % Apply Sobel Y kernel
        GyS = sum(sum(grayImage(i-1:i+1, j-1:j+1) .* SobelY));
        % Compute gradient magnitude
        SobelEdge(i,j) = sqrt(GxS^2 + GyS^2);
    end
end

% Normalize the results to display correctly
RobertEdge = uint8(RobertEdge / max(RobertEdge(:)) * 255);
PrewittEdge = uint8(PrewittEdge / max(PrewittEdge(:)) * 255);
SobelEdge = uint8(SobelEdge / max(SobelEdge(:)) * 255);

% Display the results
figure;
subplot(2,2,1); imshow(uint8(grayImage)); title('Original Image');
subplot(2,2,2); imshow(RobertEdge); title('Robert Edge Detection');
subplot(2,2,3); imshow(PrewittEdge); title('Prewitt Edge Detection');
subplot(2,2,4); imshow(SobelEdge); title('Sobel Edge Detection');


7....imagee Sharpen...

function main()
    % Load an example grayscale image
    image = imread('F:\MATLAB_DIP\labgray.jpg'); % Replace with path to your image
    if size(image, 3) == 3
        image = rgb_to_grayscale(image); % Convert to grayscale if it's an RGB image
    end

    % Apply smoothing
    smoothed_image = smooth_image(image);

    % Apply sharpening
    sharpened_image = sharpen_image(image);

    % Display original, smoothed, and sharpened images
    figure;
    subplot(1, 3, 1);
    imshow(image);
    title('Original Image');

    subplot(1, 3, 2);
    imshow(smoothed_image);
    title('Smoothed Image');

    subplot(1, 3, 3);
    imshow(sharpened_image);
    title('Sharpened Image');
end

function grayscale_image = rgb_to_grayscale(rgb_image)
    % Convert RGB image to grayscale without built-in function
    [height, width, ~] = size(rgb_image);
    grayscale_image = zeros(height, width, 'uint8');
    for i = 1:height
        for j = 1:width
            R = double(rgb_image(i, j, 1));
            G = double(rgb_image(i, j, 2));
            B = double(rgb_image(i, j, 3));
            grayscale_image(i, j) = uint8(0.2989 * R + 0.5870 * G + 0.1140 * B);
        end
    end
end

function smoothed_image = smooth_image(image)
    % Define a Gaussian kernel for smoothing (3x3)
    kernel = [1, 2, 1;
              2, 4, 2;
              1, 2, 1];
    kernel = kernel / sum(kernel(:)); % Normalize the kernel

    % Get image dimensions
    [height, width] = size(image);

    % Initialize output image
    smoothed_image = zeros(height, width, 'uint8');

    % Apply convolution without using built-in function
    for i = 2:height-1
        for j = 2:width-1
            % Extract 3x3 region
            region = double(image(i-1:i+1, j-1:j+1));
            % Apply kernel
            smoothed_value = sum(sum(region .* kernel));
            smoothed_image(i, j) = uint8(smoothed_value);
        end
    end
end

function sharpened_image = sharpen_image(image)
    % Define a sharpening kernel (Laplacian kernel for edge enhancement)
    kernel = [0, -1, 0;
             -1, 5, -1;
              0, -1, 0];

    % Get image dimensions
    [height, width] = size(image);

    % Initialize output image
    sharpened_image = zeros(height, width, 'uint8');

    % Apply convolution without using built-in function
    for i = 2:height-1
        for j = 2:width-1
            % Extract 3x3 region
            region = double(image(i-1:i+1, j-1:j+1));
            % Apply kernel
            sharpened_value = sum(sum(region .* kernel));
            % Ensure pixel value is within 0-255 range
            sharpened_image(i, j) = uint8(min(max(sharpened_value, 0), 255));
        end
    end
end

8.....Region splitting and merging......


% Read and convert the image to grayscale
inputImage = imread('input_image.jpg');  % Replace with your image file
grayImage = rgb2gray(inputImage);        % Convert to grayscale if the image is color
grayImage = double(grayImage);           % Convert to double for computation

% Parameters
threshold = 50;                          % Threshold for splitting/merging criterion (variance threshold)
minRegionSize = 8;                       % Minimum region size for merging (in pixels)
maxRegionSize = 100;                     % Maximum size to prevent splitting too much
splitThreshold = 10;                     % Variance threshold for splitting a region
regionMap = zeros(size(grayImage));      % Initialize region map

% Image dimensions
[rows, cols] = size(grayImage);

% Function to split and merge a region based on variance
function regionMap = splitMerge(image, regionMap, threshold, x, y, width, height)
    % Calculate the mean and variance of the region
    region = image(x:x+width-1, y:y+height-1);
    regionMean = mean(region(:));
    regionVar = var(region(:));

    % If variance exceeds the threshold and region is large enough, split further
    if regionVar > threshold && width > 2 && height > 2
        midX = floor(width / 2);
        midY = floor(height / 2);

        % Split region into 4 smaller sub-regions and recursively apply splitMerge
        regionMap = splitMerge(image, regionMap, threshold, x, y, midX, midY);                    % Top-left
        regionMap = splitMerge(image, regionMap, threshold, x, y + midY, midX, height - midY);     % Top-right
        regionMap = splitMerge(image, regionMap, threshold, x + midX, y, width - midX, midY);       % Bottom-left
        regionMap = splitMerge(image, regionMap, threshold, x + midX, y + midY, width - midX, height - midY);  % Bottom-right
    else
        % Mark the region as homogeneous
        regionMap(x:x+width-1, y:y+height-1) = 1;
    end
end

% Start region splitting from the full image
regionMap = splitMerge(grayImage, regionMap, splitThreshold, 1, 1, rows, cols);

% Merge small regions if their means are similar
for i = 1:rows-1
    for j = 1:cols-1
        if regionMap(i, j) == 1 && regionMap(i+1, j) == 1 && regionMap(i, j+1) == 1
            % Merge neighboring similar regions based on pixel intensity values
            if abs(grayImage(i, j) - grayImage(i+1, j)) < threshold && abs(grayImage(i, j) - grayImage(i, j+1)) < threshold
                regionMap(i, j) = 1;  % Mark as edge-free
            end
        end
    end
end

% Create the final edge-detection result
edgeImage = zeros(size(grayImage));

% Mark boundaries of regions as edges
for i = 2:rows-1
    for j = 2:cols-1
        if regionMap(i, j) ~= regionMap(i-1, j) || regionMap(i, j) ~= regionMap(i, j-1)
            edgeImage(i, j) = 255;  % Mark as boundary (edge)
        else
            edgeImage(i, j) = 0;    % Inside homogeneous regions
        end
    end
end

% Display the original image and edge-detection result
figure;
subplot(1,2,1); imshow(uint8(grayImage)); title('Original Image');
subplot(1,2,2); imshow(uint8(edgeImage)); title('Region Splitting and Merging Edge Detection');

9.....Region Growing.......

% Read the image and convert to grayscale
inputImage = imread('can1.jpg'); % Replace with your image file
grayImage = rgb2gray(inputImage);       % Convert to grayscale if the image is in color
grayImage = double(grayImage);          % Convert to double for computation

% Initialize parameters
threshold = 20;                         % Threshold for similarity
seedX = 100;                            % Seed X coordinate
seedY = 100;                            % Seed Y coordinate
[rows, cols] = size(grayImage);         % Get the size of the image
visited = false(rows, cols);            % Create a matrix to mark visited pixels
edgeImage = zeros(rows, cols);         % Initialize the edge image

% Initialize region with the seed pixel
region = zeros(rows, cols);
region(seedX, seedY) = 1;              % Mark the seed point

% Grow region starting from the seed point
while true
    changes = false;  % Flag to check if any new pixels were added

    % Loop through all pixels to expand region
    for i = 2:rows-1
        for j = 2:cols-1
            if region(i, j) == 1 && ~visited(i, j)
                visited(i, j) = true; % Mark the pixel as visited

                % Check all 8 neighbors
                for dx = -1:1
                    for dy = -1:1
                        nx = i + dx;
                        ny = j + dy;

                        % Make sure the new coordinates are within bounds
                        if nx > 0 && ny > 0 && nx <= rows && ny <= cols
                            % If the neighbor is similar to the current region, add it
                            if abs(grayImage(nx, ny) - grayImage(i, j)) < threshold && region(nx, ny) == 0
                                region(nx, ny) = 1;
                                changes = true; % New pixel added to the region
                            end
                        end
                    end
                end
            end
        end
    end

    % Stop the process if no new pixels were added in the last iteration
    if ~changes
        break;
    end
end

% Create the edge image by marking boundaries
for i = 2:rows-1
    for j = 2:cols-1
        if region(i, j) == 1
            edgeImage(i, j) = 255; % Region boundary is marked as white
        else
            edgeImage(i, j) = 0;   % Non-region areas are marked as black
        end
    end
end

% Display the results
figure;
subplot(1,2,1); imshow(uint8(grayImage)); title('Original Image');
subplot(1,2,2); imshow(uint8(edgeImage)); title('Region Growing Edge Detection');


10.....Canny....


% Read the input image
inputImage = imread('can1.jpg'); % Replace with your image file
inputImage = rgb2gray(inputImage);      % Convert to grayscale if it's a color image

% Convert the image to double for computations
inputImage = double(inputImage);

% Step 1: Apply Gaussian filter to smooth the image
% Define a 5x5 Gaussian filter
sigma = 1.4; % Standard deviation
filterSize = 5;
[x, y] = meshgrid(-floor(filterSize/2):floor(filterSize/2));
gaussianFilter = exp(-(x.^2 + y.^2) / (2 * sigma^2));
gaussianFilter = gaussianFilter / sum(gaussianFilter(:)); % Normalize

% Convolve the image with the Gaussian filter
[rows, cols] = size(inputImage);
smoothedImage = zeros(rows, cols);
for i = 3:rows-2
    for j = 3:cols-2
        region = inputImage(i-2:i+2, j-2:j+2);
        smoothedImage(i, j) = sum(sum(region .* gaussianFilter));
    end
end

% Step 2: Compute gradients using Sobel operator
Gx = [-1 0 1; -2 0 2; -1 0 1]; % Sobel x-direction kernel
Gy = [-1 -2 -1; 0 0 0; 1 2 1]; % Sobel y-direction kernel

gradientX = zeros(rows, cols);
gradientY = zeros(rows, cols);

for i = 2:rows-1
    for j = 2:cols-1
        region = smoothedImage(i-1:i+1, j-1:j+1);
        gradientX(i, j) = sum(sum(region .* Gx));
        gradientY(i, j) = sum(sum(region .* Gy));
    end
end

% Compute gradient magnitude and direction
gradientMagnitude = sqrt(gradientX.^2 + gradientY.^2);
gradientDirection = atan2d(gradientY, gradientX); % Direction in degrees

% Normalize gradient magnitude to 0-255
gradientMagnitude = gradientMagnitude / max(gradientMagnitude(:)) * 255;

% Step 3: Non-maximum suppression
nmsImage = zeros(rows, cols);
angle = mod(gradientDirection, 180); % Convert to 0-180 degrees
for i = 2:rows-1
    for j = 2:cols-1
        % Determine the neighboring pixels to compare
        if (angle(i, j) >= 0 && angle(i, j) < 22.5) || (angle(i, j) >= 157.5 && angle(i, j) <= 180)
            neighbor1 = gradientMagnitude(i, j-1);
            neighbor2 = gradientMagnitude(i, j+1);
        elseif (angle(i, j) >= 22.5 && angle(i, j) < 67.5)
            neighbor1 = gradientMagnitude(i-1, j+1);
            neighbor2 = gradientMagnitude(i+1, j-1);
        elseif (angle(i, j) >= 67.5 && angle(i, j) < 112.5)
            neighbor1 = gradientMagnitude(i-1, j);
            neighbor2 = gradientMagnitude(i+1, j);
        else
            neighbor1 = gradientMagnitude(i-1, j-1);
            neighbor2 = gradientMagnitude(i+1, j+1);
        end

        % Suppress non-maximum pixels
        if gradientMagnitude(i, j) >= neighbor1 && gradientMagnitude(i, j) >= neighbor2
            nmsImage(i, j) = gradientMagnitude(i, j);
        else
            nmsImage(i, j) = 0;
        end
    end
end

% Step 4: Double threshold and edge tracking by hysteresis
highThreshold = 40; % High threshold (adjust as needed)
lowThreshold = 20;  % Low threshold (adjust as needed)

strongEdges = nmsImage > highThreshold;
weakEdges = nmsImage > lowThreshold & nmsImage <= highThreshold;

% Edge tracking by hysteresis
finalEdges = strongEdges;
for i = 2:rows-1
    for j = 2:cols-1
        if weakEdges(i, j)
            % Check if any strong edge is in the neighborhood
            if any(any(strongEdges(i-1:i+1, j-1:j+1)))
                finalEdges(i, j) = 1;
            end
        end
    end
end

% Display the results
figure;
subplot(2, 2, 1); imshow(uint8(inputImage)); title('Original Image');
subplot(2, 2, 2); imshow(uint8(smoothedImage)); title('Smoothed Image');
subplot(2, 2, 3); imshow(uint8(nmsImage)); title('Non-Maximum Suppressed');
subplot(2, 2, 4); imshow(finalEdges); title('Final Edge Detected Image');

11....OTSu....


% Read the input image
inputImage = imread('can8.jpeg'); % Replace with your image file
inputImage = rgb2gray(inputImage);      % Convert to grayscale if it's a color image

% Convert the image to double for computations
inputImage = double(inputImage);

% Step 1: Compute Gradient Magnitude using Sobel Operator
% Define Sobel Kernels
Gx = [-1 0 1; -2 0 2; -1 0 1];
Gy = [-1 -2 -1; 0 0 0; 1 2 1];

% Initialize gradient matrices
[rows, cols] = size(inputImage);
gradientX = zeros(rows, cols);
gradientY = zeros(rows, cols);

% Apply Sobel Operator
for i = 2:rows-1
    for j = 2:cols-1
        region = inputImage(i-1:i+1, j-1:j+1);
        gradientX(i, j) = sum(sum(region .* Gx));
        gradientY(i, j) = sum(sum(region .* Gy));
    end
end

% Compute Gradient Magnitude
gradientMagnitude = sqrt(gradientX.^2 + gradientY.^2);

% Normalize gradient magnitude to 0-255
gradientMagnitude = gradientMagnitude / max(gradientMagnitude(:)) * 255;

% Step 2: Implement OTSU's Method for Thresholding
% Flatten the gradient magnitude image into a histogram
histogram = zeros(256, 1); % Initialize histogram for pixel values [0-255]

% Compute the histogram
for value = 0:255
    histogram(value + 1) = sum(gradientMagnitude(:) == value);
end

% Compute total weight and mean
totalPixels = sum(histogram);
totalMean = sum((0:255) .* histogram') / totalPixels;

% Initialize variables for OTSU
maxVariance = 0;
optimalThreshold = 0;

% Iterate through all possible thresholds
weightBackground = 0;
meanBackground = 0;

for threshold = 0:255
    weightBackground = weightBackground + histogram(threshold + 1);
    if weightBackground == 0
        continue;
    end
    weightForeground = totalPixels - weightBackground;
    if weightForeground == 0
        break;
    end

    meanBackground = meanBackground + threshold * histogram(threshold + 1);
    currentMeanBackground = meanBackground / weightBackground;
    currentMeanForeground = (totalMean - meanBackground) / weightForeground;

    % Compute Between-Class Variance
    betweenClassVariance = weightBackground * weightForeground * ...
        (currentMeanBackground - currentMeanForeground)^2;

    % Update maximum variance and optimal threshold
    if betweenClassVariance > maxVariance
        maxVariance = betweenClassVariance;
        optimalThreshold = threshold;
    end
end

% Step 3: Threshold the Gradient Magnitude
binaryEdgeImage = gradientMagnitude > optimalThreshold;

% Display the Results
figure;
subplot(1, 3, 1); imshow(uint8(inputImage)); title('Original Image');
subplot(1, 3, 2); imshow(uint8(gradientMagnitude)); title('Gradient Magnitude');
subplot(1, 3, 3); imshow(binaryEdgeImage); title('Edge Detected Image (OTSU)');

12....Background subtraction......

% Background Subtraction in MATLAB

% Read the background image
background = imread('otsu.jpg');
background = rgb2gray(background); % Convert to grayscale if it's RGB

% Read the current frame (e.g., from a video or an image)
current_frame = imread('graynegative.jpg');
current_frame = rgb2gray(current_frame); % Convert to grayscale if it's RGB

% Convert images to double precision for subtraction
background = double(background);
current_frame = double(current_frame);

% Perform background subtraction
difference = abs(current_frame - background);

% Thresholding the difference to identify foreground objects
threshold = 30; % Adjust based on your needs
foreground = difference > threshold;

% Display the original, background, and background-subtracted images
subplot(1, 3, 1);
imshow(uint8(background));
title('Background');

subplot(1, 3, 2);
imshow(uint8(current_frame));
title('Current Frame');

subplot(1, 3, 3);
imshow(foreground);
title('Background Subtracted Image');

% Optional: Save the result
imwrite(foreground, 'background_subtracted_image.jpg');

......filter noise enhance....

% Load an image
originalImage = imread('pepper.jpg'); % Load your image file
if size(originalImage, 3) == 3
    originalImage = rgb2gray(originalImage); % Convert to grayscale for simplicity
end

% Resize the image for faster processing
originalImage = imresize(originalImage, [256, 256]);

% Display original image
figure;
subplot(3, 3, 1);
imshow(originalImage);
title('Original Image');

% Apply different types of noise
noisyImages = {};
noiseTypes = {'gaussian', 'salt & pepper', 'speckle'};

for i = 1:length(noiseTypes)
    noisyImages{i} = imnoise(originalImage, noiseTypes{i});
    subplot(3, 3, i+1);
    imshow(noisyImages{i});
    title([noiseTypes{i}, ' Noise']);
end

% Restore/enhance the noisy images using different filters
restoredImages = {};
filterMethods = {'Gaussian', 'Median', 'Wiener'};
PSNRValues = zeros(length(noiseTypes), length(filterMethods)); % For performance comparison

for i = 1:length(noiseTypes)
    for j = 1:length(filterMethods)
        switch filterMethods{j}
            case 'Gaussian'
                h = fspecial('gaussian', [3, 3], 0.5);
                restoredImages{i, j} = imfilter(noisyImages{i}, h, 'replicate');
            case 'Median'
                restoredImages{i, j} = medfilt2(noisyImages{i}, [3, 3]);
            case 'Wiener'
                restoredImages{i, j} = wiener2(noisyImages{i}, [5, 5]);
        end
        % Compute PSNR for performance evaluation
        PSNRValues(i, j) = psnr(restoredImages{i, j}, originalImage);
    end
end

% Display restored images and PSNR values
figure;
for i = 1:length(noiseTypes)
    for j = 1:length(filterMethods)
        subplot(length(noiseTypes), length(filterMethods), (i-1)*length(filterMethods) + j);
        imshow(restoredImages{i, j});
        title(sprintf('%s Filter\n(PSNR: %.2f dB)', filterMethods{j}, PSNRValues(i, j)));
    end
end

% Print PSNR values
disp('PSNR Values (dB):');
disp(array2table(PSNRValues, 'VariableNames', filterMethods, 'RowNames', noiseTypes));