Calculating UACI, NPCR for 2 images

Question

LOKESH 2012-4-22

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/36171-calculating-uaci-npcr-for-2-images

评论： sajjad 2025-3-16

Presently I have a problem in implementing the UACI and PSNR code for 2 images which are respectively Original & Cipher Images. I need the code for calculating it......

What would the Code? I don't have much idea about it.. Do I need to find UACI between 2 images or only 1 image will be sufficient

5 个评论
显示 3更早的评论隐藏 3更早的评论

sajjad 2025-3-16

Professor, i am using a choatic map (in code it is with name a squence, this sequence is gotten from a matrix S2, that S2 (minima of a chaotic system, i have already solved that chatic system no need to solve it again) ) for image encrotion by using algorithms given in this paper (Image encryption using chaos and block cipher, with DOI No10.5539/cis.v4n1p172). But code is not running corecctly. Please checkt it. I have attached the all files. It is an image encryption scheme based on combination of pixel shuffling and new modified version of simplified AES. Chaotic map is used for shuffling and improving S-AES efficiency through S-box design. Chaos is used to expand diffusion and confusion in the image.

'' in the paper he use a baker's map "

(FOLLOWIGN IS THE CODE FOR MY MAP i,e SEQUENCE instead OF BAKER's MAP)

%% (1)

% Function for Shuffling Pixels Using Chaotic Sequence (Replacing Baker's Map)

function shuffled_img = shuffle_pixels_with_sequence(img, sequence)

[rows, cols] = size(img);

num_points = size(sequence, 1);

% Generate initial grid coordinates

[X, Y] = meshgrid(1:cols, 1:rows);

X = X(:) / cols;

Y = Y(:) / rows;

% Iterate to compute the shuffled positions using the chaotic sequence

for i = 1:num_points

X_next = mod(X + sequence(i, 1), 1);

Y_next = mod(Y + sequence(i, 2), 1);

X = X_next;

Y = Y_next;

end

% Ensure X and Y values are within the valid range

X = ceil(X * cols);

Y = ceil(Y * rows);

% Make sure X and Y are within the range [1, rows] and [1, cols]

X = max(min(X, cols), 1);

Y = max(min(Y, rows), 1);

% Rearrange pixels based on shuffled indices

indices = sub2ind([rows, cols], Y, X);

shuffled_img = img(indices);

shuffled_img = reshape(shuffled_img, rows, cols);

% **Step 1: Vertical Permutation**

Pmap = randperm(rows);

shuffled_img = shuffled_img(Pmap, :);

% **Step 2: Horizontal Permutation**

Pmap = randperm(cols);

shuffled_img = shuffled_img(:, Pmap);

% **Step 3: Apply Circular Shift for Diagonal Permutation**

for it = 2:rows

Shiftsize = mod(num_points - it + 1, cols); % Prevents exceeding bounds

shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, Shiftsize]);

end

% **Step 4: Apply Vertical Permutation Again**

shuffled_img = shuffled_img(Pmap, :);

% **Step 5: Apply Final Diagonal Permutation**

for it = 2:rows

Shiftsize = mod(num_points - it + 1, cols);

shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, -Shiftsize]);

end

%% (2)

function encrypted_img = s_aes_encrypt(img, sbox)

[rows, cols] = size(img);

% Ensure the S-box is expanded to match the image size

if size(sbox, 1) < rows || size(sbox, 2) < cols

sbox_expanded = repmat(sbox, ceil(rows / 256), ceil(cols / 256));

sbox_expanded = sbox_expanded(1:rows, 1:cols);

else

sbox_expanded = sbox(1:rows, 1:cols);

end

% Convert image to double for computations

img = double(img);

% Round 1

img = add_round_key(img, sbox_expanded); % Add Round Key

img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles

img = shift_rows(img); % Shift Rows

img = mix_columns(img); % Mix Columns

img = add_round_key(img, sbox_expanded); % Add Round Key

% Round 2

img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles

img = shift_rows(img); % Shift Rows

img = add_round_key(img, sbox_expanded); % Add Round Key

% Convert back to uint8 for image display

encrypted_img = uint8(img);

end

%% (3)

function state = substitute_nibbles(state, sbox)

[rows, cols] = size(state);

% Substitute each nibble using the S-box

for i = 1:rows

for j = 1:cols

% Ensure the value is within the range of the S-box

nibble_value = mod(state(i, j), 256) + 1; % Map to 1-256

state(i, j) = sbox(nibble_value);

end

%% (4)

function state = shift_rows(state)

% Shift rows as per S-AES

state(2, :) = circshift(state(2, :), -1); % Shift second row by 1 position

end

%%% (5)

function state = mix_columns(state)

% Mix columns as per S-AES

temp = state;

state(1, 1) = bitxor(temp(1, 1), temp(2, 1));

state(1, 2) = bitxor(temp(1, 2), temp(2, 2));

state(2, 1) = bitxor(temp(1, 1), temp(2, 1));

state(2, 2) = bitxor(temp(1, 2), temp(2, 2));

end

%%% (6)

function state = add_round_key(state, round_key)

% XOR the state with the round key

state = bitxor(state, round_key);

end

%%%%(7)

%% (5) generate_sbox - Fully Corrected Version

function sbox = generate_sbox(sequence)

NoIt = length(sequence);

D = zeros(1, NoIt);

% **Step 1-5: Compute D(it) Using Iterative Transformation**

for it = 1:NoIt

idx = mod(it - 1, length(sequence)) + 1;

X = sequence(idx, 1);

Y = sequence(idx, 2);

% ? Corrected: Use iterative transformations

X_next = mod(X + Y, 1);

Y_next = mod(X + 2 * Y, 1);

% Update X, Y for next iteration

X = X_next;

Y = Y_next;

% ? Compute D(it) using correct transformation

D(it) = sqrt(abs(X.^3 + Y.^3));

end

% **Step 6-10: Apply Sorting & Nonlinear Transformation**

D = sort(D, 'ascend');

M = max(D);

for it1 = 1:NoIt

for it2 = 1:NoIt

if D(it2) == M

D(it2) = -16 + it1;

end

% **Step 11-20: Normalize & Create the S-Box**

D = abs(D);

max_D = max(D);

sbox = uint8(mod(round(D * 255 / max_D), 256));

% **Ensure the S-Box is exactly (256 × 256)**

sbox = repmat(sbox(:), ceil(256*256/numel(sbox)), 1);

sbox = reshape(sbox(1:256*256), 256, 256);

end

%%%%(8)

%% (8) encryption_main_subscript - Fully Matches Figure 5

clc; clear; close all;

%% **Step 1: Load and Preprocess Image**

filePath = 'C:\Users\HP EliteBook 840 G3\Desktop\pic3.jpeg';

img = imread(filePath);

% Convert to grayscale if the image is RGB

if size(img, 3) == 3

img = rgb2gray(img);

end

% Resize image to 256x256 (if not already)

img = imresize(img, [256, 256]);

% Display (a) Original Image

figure;

subplot(2, 3, 1);

imshow(img);

title('(a) Original Image');

%% **Step 2: Load Sequence from S2.xlsx for Pixel Shuffling**

sequenceFile = 'C:\Users\HP EliteBook 840 G3\Desktop\S1.xlsx';

data = xlsread(sequenceFile);

data = data(:); % Flatten matrix

% Scale the sequence between 0 and 1 (Normalization)

scaled_sequence = (data - min(data)) / (max(data) - min(data));

maxima_x = scaled_sequence(1:end-1);

maxima_y = scaled_sequence(2:end);

sequence = [maxima_x, maxima_y];

%% **Step 3: Apply Algorithm 1 (Pixel Shuffling)**

shuffledImg = shuffle_pixels_with_sequence(img, sequence);

% Display (b) Shuffled Image

subplot(2, 3, 2);

imshow(shuffledImg);

title('(b) Shuffled Image');

%% **Step 4: Apply Algorithm 2 (Generate Key-Dependent S-Box)**

sbox = generate_sbox(sequence); % Generate dynamic S-Box

%% **Step 5: Apply Simplified AES (S-AES) Encryption**

I = s_aes_encrypt(shuffledImg, sbox); % Pass S-Box explicitly

% Display (c) Ciphered Image

subplot(2, 3, 3);

imshow(I);

title('(c) Encrypted Image');

%% **Step 6: Compute and Display Histograms**

subplot(2, 3, 4); imhist(img); title('(d) Histogram of Original Image');

subplot(2, 3, 5); imhist(shuffledImg); title('(e) Histogram of Shuffled Image');

subplot(2, 3, 6); imhist(I); title('(f) Histogram of Encrypted Image');

% %% **Step 7: Save Encrypted Image**

% outputPath = 'C:\Users\HP EliteBook 840 G3\Desktop\encrypted_image.png';

% imwrite(cipheredImg, outputPath);

% disp(['Encrypted image saved to: ', outputPath]);

sajjad 2025-3-16

For S1; from the following code just consider MAXIMA and in above codes i have converted S1 in sequence.

The system this paper (DOI: 10.1017/S0022112003003835) Thank you.

%% funn

function dydt = funn(t,y,x,nx,dx,delta,Reynolds,H,Hdot);

%persistent iter

%if isempty(iter)

% iter = 0;

%end

%iter = iter + 1;

%if mod(iter,1000)==0

% t

% iter = 0;

%end

G = y(1:nx);

F = y(nx+1:2*nx);

% F

% Compute spatial derivatives

dFdx = zeros(nx,1);

d2Fdx2 = zeros(nx,1);

dFdx(2:nx-1) = (F(3:nx)-F(1:nx-2))/(2*dx);

dFdx(nx) = Hdot(t); % This corresponds to F'(1,t) = H'(t)

d2Fdx2(2:nx-1) = (F(3:nx)-2*F(2:nx-1)+F(1:nx-2))/dx^2;

%d2Fdx2(nx) = (dFdx(nx)-dFdx(nx-1))/dx;

Fnxp1 = F(nx-1) + 2*dx*dFdx(nx);

d2Fdx2(nx) = (Fnxp1-2*F(nx)+F(nx-1))/dx^2;

% Compute temporal derivatives

dFdt = zeros(nx,1);

dFdt(1) = F(1);

dFdt(2:nx-1) = d2Fdx2(2:nx-1)+dFdx(2:nx-1)./x(2:nx-1)-...

F(2:nx-1)./x(2:nx-1).^2+H(t)^2*G(2:nx-1);

dFdt(nx) = F(nx) + Hdot(t); % This corresponds to F(1,t) =- H'(t)

% G

% Compute spatial derivatives

dGdx = zeros(nx,1);

d2Gdx2 = zeros(nx,1);

dGdx(2:nx-1) = (G(3:nx)-G(1:nx-2))/(2*dx);

d2Gdx2(2:nx-1) = (G(3:nx)-2*G(2:nx-1)+G(1:nx-2))/dx^2;

% Compute temporal derivatives

dGdt = zeros(nx,1);

dGdt(1) = G(1);

dGdt(2:nx-1) = Hdot(t)/H(t).*x(2:nx-1).*dGdx(2:nx-1)+ ...

1/H(t).*F(2:nx-1).*dGdx(2:nx-1)- ...

1/H(t).*dFdx(2:nx-1).*G(2:nx-1)- ...

2.*F(2:nx-1).*G(2:nx-1)./(H(t)*x(2:nx-1))+ ...

(d2Gdx2(2:nx-1)+dGdx(2:nx-1)./x(2:nx-1)-G(2:nx-1)./x(2:nx-1).^2)./ ...

(H(t)^2*Reynolds);

dGdt(nx) = G(nx) - (d2Fdx2(nx)+dFdx(nx)/x(nx)-F(nx)/x(nx)^2)/(-H(t)^2);

%Taken from the article, should be the same as set

%dGdt(nx) = G(nx) + 2/H(t)^2*((F(nx-1)+Hdot(t)*(1+dx))/dx^2 + Hdot(t));

dydt = [dGdt;dFdt];

end

%%% MAIN SUBSCRIPTT

format long;

clear all;

close all;

clc;

% Parameters

delta = 0.3; % Amplitude of H

Reynolds = 900; % Reynolds number

tstart = 0; % Start time

tend = 206400; % End time

nx = 101; % Number of spatial points

xstart = 0.0; % Spatial domain start

xend = 1.0; % Spatial domain end

x = linspace(xstart, xend, nx).'; % Spatial grid

dx = (xend - xstart) / (nx - 1); % Spatial step size

tspan = [tstart, tend]; % Time span

% numpoints = 1000;

% tspan = linspace(tstart, tend,numpoints); % Time span

% Initial conditions

G0 = zeros(nx, 1); % Initial condition for G

F0 = zeros(nx, 1); % Initial condition for F

y0 = [G0; F0]; % Combine into one vector

% Define H(t) and its derivative

H = @(t) 1 + delta * cos(2 * t);

Hdot = @(t) -2 * delta * sin(2 * t);

% Mass matrix for DAE

M = [eye(nx), zeros(nx); zeros(nx, 2*nx)];

M(1, 1) = 0; % Boundary condition G(0,t) = 0

M(nx, nx) = 0; % Boundary condition G(1,t) = 0

options = odeset('Mass', M, 'RelTol', 1e-3, 'AbsTol', 1e-6);

% Solve the system

[T, Y] = ode23t(@(t, y) funn(t, y, x, nx, dx, delta, Reynolds, H, Hdot), tspan, y0, options);

% Extract G and F

G = Y(:, 1:nx);

F = Y(:, nx+1:2*nx);

% Evaluate G at eta = 1 using interpolation

eta_query = 1; % Change from 1/4 to 1

G_eta_1 = zeros(size(T)); % Preallocate

for i = 1:length(T)

G_eta_1(i) = interp1(x, G(i, :), eta_query, 'spline'); % Interpolation

end

% Plot G(1, t) over time

figure(1);

plot(T, G_eta_1, 'LineWidth', 1.5);

xlabel('Time (t)');

ylabel('G(1, t)');

title('Time evolution of G(1, t)');

grid on;

[maxima_G_eta_1, locs_max_G_eta_1] = findpeaks(G_eta_1); % Maxima of G(1/4, t)

[minima_G_eta_1, locs_min_G_eta_1] = findpeaks(-G_eta_1); % Minima of G(1/4, t)

minima_G_eta_1 = -minima_G_eta_1; % Correct sign of minima

% Display counts of maxima and minima

sajjad1 = length(maxima_G_eta_1);

sajjad2 = length(minima_G_eta_1);

disp(['Number of maxima: ', num2str(sajjad1)]);

disp(['Number of minima: ', num2str(sajjad2)]);

% Multiply maxima and minima by 10^15 and apply mod 256

scaled_maxima = mod(round(maxima_G_eta_1 * 1e15), 256);

scaled_minima = mod(round(minima_G_eta_1 * 1e15), 256);

% Display scaled maxima and minima

disp('Scaled maxima (mod 256):');

disp(scaled_maxima);

disp('Scaled minima (mod 256):');

disp(scaled_minima);

% Return Map for Maxima

figure(2); % New figure

plot(maxima_G_eta_1(1:end-1), maxima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'r', 'MarkerSize', 4);

xlabel('Maxima_{n}');

ylabel('Maxima_{n+1}');

title('Return Map of Maxima');

grid on;

% Return Map for Minima

figure(3); % Another new figure

plot(minima_G_eta_1(1:end-1), minima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'g', 'MarkerSize', 4);

xlabel('Minima_{n}');

ylabel('Minima_{n+1}');

title('Return Map of Minima');

grid on;

Please feel to ask any other information about it.

Thank you.

sajjad 2025-3-16

Professor, i am using a choatic map (in code it is with name a squence, this sequence is gotten from a matrix S2, that S2 (minima of a chaotic system, i have already solved that chatic system no need to solve it again) ) for image encrotion by using algorithms given in this paper (Image encryption using chaos and block cipher, with DOI No10.5539/cis.v4n1p172). But code is not running corecctly. Please checkt it. I have attached the all files. It is an image encryption scheme based on combination of pixel shuffling and new modified version of simplified AES. Chaotic map is used for shuffling and improving S-AES efficiency through S-box design. Chaos is used to expand diffusion and confusion in the image.

'' in the paper he use a baker's map "

(FOLLOWIGN IS THE CODE FOR MY MAP i,e SEQUENCE instead OF BAKER's MAP)

%% (1)

% Function for Shuffling Pixels Using Chaotic Sequence (Replacing Baker's Map by my sequence)

function shuffled_img = shuffle_pixels_with_sequence(img, sequence)

[rows, cols] = size(img);

num_points = size(sequence, 1);

% Generate initial grid coordinates

[X, Y] = meshgrid(1:cols, 1:rows);

X = X(:) / cols;

Y = Y(:) / rows;

% Iterate to compute the shuffled positions using the chaotic sequence

for i = 1:num_points

X_next = mod(X + sequence(i, 1), 1);

Y_next = mod(Y + sequence(i, 2), 1);

X = X_next;

Y = Y_next;

end

% Ensure X and Y values are within the valid range

X = ceil(X * cols);

Y = ceil(Y * rows);

% Make sure X and Y are within the range [1, rows] and [1, cols]

X = max(min(X, cols), 1);

Y = max(min(Y, rows), 1);

% Rearrange pixels based on shuffled indices

indices = sub2ind([rows, cols], Y, X);

shuffled_img = img(indices);

shuffled_img = reshape(shuffled_img, rows, cols);

% **Step 1: Vertical Permutation**

Pmap = randperm(rows);

shuffled_img = shuffled_img(Pmap, :);

% **Step 2: Horizontal Permutation**

Pmap = randperm(cols);

shuffled_img = shuffled_img(:, Pmap);

% **Step 3: Apply Circular Shift for Diagonal Permutation**

for it = 2:rows

Shiftsize = mod(num_points - it + 1, cols); % Prevents exceeding bounds

shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, Shiftsize]);

end

% **Step 4: Apply Vertical Permutation Again**

shuffled_img = shuffled_img(Pmap, :);

% **Step 5: Apply Final Diagonal Permutation**

for it = 2:rows

Shiftsize = mod(num_points - it + 1, cols);

shuffled_img(it, :) = circshift(shuffled_img(it, :), [0, -Shiftsize]);

end

%% (2)

function encrypted_img = s_aes_encrypt(img, sbox)

[rows, cols] = size(img);

% Ensure the S-box is expanded to match the image size

if size(sbox, 1) < rows || size(sbox, 2) < cols

sbox_expanded = repmat(sbox, ceil(rows / 256), ceil(cols / 256));

sbox_expanded = sbox_expanded(1:rows, 1:cols);

else

sbox_expanded = sbox(1:rows, 1:cols);

end

% Convert image to double for computations

img = double(img);

% Round 1

img = add_round_key(img, sbox_expanded); % Add Round Key

img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles

img = shift_rows(img); % Shift Rows

img = mix_columns(img); % Mix Columns

img = add_round_key(img, sbox_expanded); % Add Round Key

% Round 2

img = substitute_nibbles(img, sbox_expanded); % Substitute Nibbles

img = shift_rows(img); % Shift Rows

img = add_round_key(img, sbox_expanded); % Add Round Key

% Convert back to uint8 for image display

encrypted_img = uint8(img);

end

%% (3)

function state = substitute_nibbles(state, sbox)

[rows, cols] = size(state);

% Substitute each nibble using the S-box

for i = 1:rows

for j = 1:cols

% Ensure the value is within the range of the S-box

nibble_value = mod(state(i, j), 256) + 1; % Map to 1-256

state(i, j) = sbox(nibble_value);

end

%% (4)

function state = shift_rows(state)

% Shift rows as per S-AES

state(2, :) = circshift(state(2, :), -1); % Shift second row by 1 position

end

%% (5)

function state = mix_columns(state)

% Mix columns as per S-AES

temp = state;

state(1, 1) = bitxor(temp(1, 1), temp(2, 1));

state(1, 2) = bitxor(temp(1, 2), temp(2, 2));

state(2, 1) = bitxor(temp(1, 1), temp(2, 1));

state(2, 2) = bitxor(temp(1, 2), temp(2, 2));

end

%% (6)

function state = add_round_key(state, round_key)

% XOR the state with the round key

state = bitxor(state, round_key);

end

%% (7)

%% (5) generate_sbox - Fully Corrected Version

function sbox = generate_sbox(sequence)

NoIt = length(sequence);

D = zeros(1, NoIt);

% **Step 1-5: Compute D(it) Using Iterative Transformation**

for it = 1:NoIt

idx = mod(it - 1, length(sequence)) + 1;

X = sequence(idx, 1);

Y = sequence(idx, 2);

% ? Corrected: Use iterative transformations

X_next = mod(X + Y, 1);

Y_next = mod(X + 2 * Y, 1);

% Update X, Y for next iteration

X = X_next;

Y = Y_next;

% ? Compute D(it) using correct transformation

D(it) = sqrt(abs(X.^3 + Y.^3));

end

% **Step 6-10: Apply Sorting & Nonlinear Transformation**

D = sort(D, 'ascend');

M = max(D);

for it1 = 1:NoIt

for it2 = 1:NoIt

if D(it2) == M

D(it2) = -16 + it1;

end

% **Step 11-20: Normalize & Create the S-Box**

D = abs(D);

max_D = max(D);

sbox = uint8(mod(round(D * 255 / max_D), 256));

% **Ensure the S-Box is exactly (256 × 256)**

sbox = repmat(sbox(:), ceil(256*256/numel(sbox)), 1);

sbox = reshape(sbox(1:256*256), 256, 256);

end

%%% (8)

%% (8) encryption_main_subscript - Fully Matches Figure 5

clc; clear; close all;

%% **Step 1: Load and Preprocess Image**

filePath = 'C:\Users\HP EliteBook 840 G3\Desktop\pic1.jpeg';

img = imread(filePath);

% Convert to grayscale if the image is RGB

if size(img, 3) == 3

img = rgb2gray(img);

end

% Resize image to 256x256 (if not already)

img = imresize(img, [256, 256]);

% Display (a) Original Image

figure;

subplot(2, 3, 1);

imshow(img);

title('(a) Original Image');

%% **Step 2: Load Sequence from S2.xlsx for Pixel Shuffling**

sequenceFile = 'C:\Users\HP EliteBook 840 G3\Desktop\S2.xlsx';

data = xlsread(sequenceFile);

data = data(:); % Flatten matrix

% Scale the sequence between 0 and 1 (Normalization)

scaled_sequence = (data - min(data)) / (max(data) - min(data));

maxima_x = scaled_sequence(1:end-1);

maxima_y = scaled_sequence(2:end);

sequence = [maxima_x, maxima_y];

%% **Step 3: Apply Algorithm 1 (Pixel Shuffling)**

shuffledImg = shuffle_pixels_with_sequence(img, sequence);

% Display (b) Shuffled Image

subplot(2, 3, 2);

imshow(shuffledImg);

title('(b) Shuffled Image');

%% **Step 4: Apply Algorithm 2 (Generate Key-Dependent S-Box)**

sbox = generate_sbox(sequence); % Generate dynamic S-Box

%% **Step 5: Apply Simplified AES (S-AES) Encryption**

I = s_aes_encrypt(shuffledImg, sbox); % Pass S-Box explicitly

% Display (c) Ciphered Image

subplot(2, 3, 3);

imshow(I);

title('(c) Encrypted Image');

%% **Step 6: Compute and Display Histograms**

subplot(2, 3, 4); imhist(img); title('(d) Histogram of Original Image');

subplot(2, 3, 5); imhist(shuffledImg); title('(e) Histogram of Shuffled Image');

subplot(2, 3, 6); imhist(I); title('(f) Histogram of Encrypted Image');

% %% **Step 7: Save Encrypted Image**

% outputPath = 'C:\Users\HP EliteBook 840 G3\Desktop\encrypted_image.png';

% imwrite(cipheredImg, outputPath);

% disp(['Encrypted image saved to: ', outputPath]);

sajjad 2025-3-16

For S1; from the following code just consider MAXIMA and in above codes i have converted S1 in sequence.

The system this paper (DOI: 10.1017/S0022112003003835) Thank you.

%% funn

function dydt = funn(t,y,x,nx,dx,delta,Reynolds,H,Hdot);

%persistent iter

%if isempty(iter)

% iter = 0;

%end

%iter = iter + 1;

%if mod(iter,1000)==0

% t

% iter = 0;

%end

G = y(1:nx);

F = y(nx+1:2*nx);

% F

% Compute spatial derivatives

dFdx = zeros(nx,1);

d2Fdx2 = zeros(nx,1);

dFdx(2:nx-1) = (F(3:nx)-F(1:nx-2))/(2*dx);

dFdx(nx) = Hdot(t); % This corresponds to F'(1,t) = H'(t)

d2Fdx2(2:nx-1) = (F(3:nx)-2*F(2:nx-1)+F(1:nx-2))/dx^2;

%d2Fdx2(nx) = (dFdx(nx)-dFdx(nx-1))/dx;

Fnxp1 = F(nx-1) + 2*dx*dFdx(nx);

d2Fdx2(nx) = (Fnxp1-2*F(nx)+F(nx-1))/dx^2;

% Compute temporal derivatives

dFdt = zeros(nx,1);

dFdt(1) = F(1);

dFdt(2:nx-1) = d2Fdx2(2:nx-1)+dFdx(2:nx-1)./x(2:nx-1)-...

F(2:nx-1)./x(2:nx-1).^2+H(t)^2*G(2:nx-1);

dFdt(nx) = F(nx) + Hdot(t); % This corresponds to F(1,t) =- H'(t)

% G

% Compute spatial derivatives

dGdx = zeros(nx,1);

d2Gdx2 = zeros(nx,1);

dGdx(2:nx-1) = (G(3:nx)-G(1:nx-2))/(2*dx);

d2Gdx2(2:nx-1) = (G(3:nx)-2*G(2:nx-1)+G(1:nx-2))/dx^2;

% Compute temporal derivatives

dGdt = zeros(nx,1);

dGdt(1) = G(1);

dGdt(2:nx-1) = Hdot(t)/H(t).*x(2:nx-1).*dGdx(2:nx-1)+ ...

1/H(t).*F(2:nx-1).*dGdx(2:nx-1)- ...

1/H(t).*dFdx(2:nx-1).*G(2:nx-1)- ...

2.*F(2:nx-1).*G(2:nx-1)./(H(t)*x(2:nx-1))+ ...

(d2Gdx2(2:nx-1)+dGdx(2:nx-1)./x(2:nx-1)-G(2:nx-1)./x(2:nx-1).^2)./ ...

(H(t)^2*Reynolds);

dGdt(nx) = G(nx) - (d2Fdx2(nx)+dFdx(nx)/x(nx)-F(nx)/x(nx)^2)/(-H(t)^2);

%Taken from the article, should be the same as set

%dGdt(nx) = G(nx) + 2/H(t)^2*((F(nx-1)+Hdot(t)*(1+dx))/dx^2 + Hdot(t));

dydt = [dGdt;dFdt];

end

%%% MAIN SUBSCRIPTT

format long;

clear all;

close all;

clc;

% Parameters

delta = 0.3; % Amplitude of H

Reynolds = 900; % Reynolds number

tstart = 0; % Start time

tend = 206400; % End time

nx = 101; % Number of spatial points

xstart = 0.0; % Spatial domain start

xend = 1.0; % Spatial domain end

x = linspace(xstart, xend, nx).'; % Spatial grid

dx = (xend - xstart) / (nx - 1); % Spatial step size

tspan = [tstart, tend]; % Time span

% numpoints = 1000;

% tspan = linspace(tstart, tend,numpoints); % Time span

% Initial conditions

G0 = zeros(nx, 1); % Initial condition for G

F0 = zeros(nx, 1); % Initial condition for F

y0 = [G0; F0]; % Combine into one vector

% Define H(t) and its derivative

H = @(t) 1 + delta * cos(2 * t);

Hdot = @(t) -2 * delta * sin(2 * t);

% Mass matrix for DAE

M = [eye(nx), zeros(nx); zeros(nx, 2*nx)];

M(1, 1) = 0; % Boundary condition G(0,t) = 0

M(nx, nx) = 0; % Boundary condition G(1,t) = 0

options = odeset('Mass', M, 'RelTol', 1e-3, 'AbsTol', 1e-6);

% Solve the system

[T, Y] = ode23t(@(t, y) funn(t, y, x, nx, dx, delta, Reynolds, H, Hdot), tspan, y0, options);

% Extract G and F

G = Y(:, 1:nx);

F = Y(:, nx+1:2*nx);

% Evaluate G at eta = 1 using interpolation

eta_query = 1; % Change from 1/4 to 1

G_eta_1 = zeros(size(T)); % Preallocate

for i = 1:length(T)

G_eta_1(i) = interp1(x, G(i, :), eta_query, 'spline'); % Interpolation

end

% Plot G(1, t) over time

figure(1);

plot(T, G_eta_1, 'LineWidth', 1.5);

xlabel('Time (t)');

ylabel('G(1, t)');

title('Time evolution of G(1, t)');

grid on;

[maxima_G_eta_1, locs_max_G_eta_1] = findpeaks(G_eta_1); % Maxima of G(1/4, t)

[minima_G_eta_1, locs_min_G_eta_1] = findpeaks(-G_eta_1); % Minima of G(1/4, t)

minima_G_eta_1 = -minima_G_eta_1; % Correct sign of minima

% Display counts of maxima and minima

sajjad1 = length(maxima_G_eta_1);

sajjad2 = length(minima_G_eta_1);

disp(['Number of maxima: ', num2str(sajjad1)]);

disp(['Number of minima: ', num2str(sajjad2)]);

% Multiply maxima and minima by 10^15 and apply mod 256

scaled_maxima = mod(round(maxima_G_eta_1 * 1e15), 256);

scaled_minima = mod(round(minima_G_eta_1 * 1e15), 256);

% Display scaled maxima and minima

disp('Scaled maxima (mod 256):');

disp(scaled_maxima);

disp('Scaled minima (mod 256):');

disp(scaled_minima);

% Return Map for Maxima

figure(2); % New figure

plot(maxima_G_eta_1(1:end-1), maxima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'r', 'MarkerSize', 4);

xlabel('Maxima_{n}');

ylabel('Maxima_{n+1}');

title('Return Map of Maxima');

grid on;

% Return Map for Minima

figure(3); % Another new figure

plot(minima_G_eta_1(1:end-1), minima_G_eta_1(2:end), 'o','LineWidth', 2,'Color', 'g', 'MarkerSize', 4);

xlabel('Minima_{n}');

ylabel('Minima_{n+1}');

title('Return Map of Minima');

grid on;

Please feel to ask any other information about it.

Thank you.

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Image Analyst 2012-4-22

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/36171-calculating-uaci-npcr-for-2-images#answer_45266

在 MATLAB Online 中打开

See my demo:

% Demo to calculate PSNR of a gray scale image.
% http://en.wikipedia.org/wiki/PSNR
% by ImageAnalyst 
clc; 
close all; 
clearvars; 
workspace;
% Read in standard MATLAB demo image.
grayImage = imread('cameraman.tif');
[rows columns] = size(grayImage);
subplot(2, 2, 1);
imshow(grayImage, []);
title('Original Grey Scale Image');
set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.
% Add noise to it.
noisyImage = imnoise(grayImage, 'gaussian', 0, 0.003);
subplot(2, 2, 2);
imshow(noisyImage, []);
title('Noisy Image');
% Calculate mean square error.
mseImage = (double(grayImage) - double(noisyImage)) .^ 2;
subplot(2, 2, 3);
imshow(mseImage, []);
title('MSE Image');
mse = sum(sum(mseImage)) / (rows * columns);
% Calculate PSNR (Peak Signal to noise ratio).
PSNR = 10 * log10( 256^2 / mse);
message = sprintf('The mean square error is %.2f\nThe PSNR = %.2f', ...
  mse, PSNR);
msgbox(message);
% set(gcf, 'Position', get(0,'Screensize')); % Maximize figure.

After this example you should be able to use similar code to easily calculate NPCR and UACI, according to http://www.cyberjournals.com/Papers/Apr2011/05.pdf. You need both images to calculate them since they are a measure of the differences between the two images.

2 个评论
显示无隐藏无

Komeil 2016-10-20

The NPCR is different from the PSNR.

Walter Roberson 2016-10-20

Image Analyst says "After this example you should be able to use similar code to easily calculate NPCR" -- in other words, he has not given the code here, but he has given a framework that can be adapted to calculate NPCR by anyone who looks up the formula for NPCR.

请先登录，再进行评论。

Answer 2

Hyderkkk 2020-8-31

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/36171-calculating-uaci-npcr-for-2-images#answer_487220

How do i get the NPCR and UACI of an original speech and encryption in matlab code ?

14 个评论
显示 12更早的评论隐藏 12更早的评论

Walter Roberson 2021-5-25

This file give the relation between two encrypted images of same size and type but I need the relation between original speech and encrypted speech

Are you familiar with the encryption technique that is often called ROT13? It is a caeser cypher in which each letter in the first half of the alphabet is replaced with the letter in the corresponding position in the second half of the alphabet, and vice versa, so 'A' gets replaced with the character 13 further on (which is 'N'), and 'p' would get replaced with the character 13 earlier (which is 'c')

Now take the original image, and apply ROT13 to it. Then take the result, and apply ROT13 to that, for even stronger encryption.

Now take the relationship between that double-ROT13 encrypted image and the encryption you found using a different method. As they are both encrypted image, there should be no problem comparing the two using an algorithm that compares encrypted images, right?

... But it turns out that applying ROT13 twice always gets you back exactly the original data. Therefore the original data is itself an encrypted version of the original data... just encrypted with a very weak encryption. And so you can use your algorithm that compares two encrypted images in order to compare the original image to the version of it you encrypted some other way than double-ROT13.

sajjad 2025-1-24

I have almost completed the image encryption. Now how do i calculate the NPCR and UACI of an original image and encryption in matlab code ?

Please guide me in this regard.

Thank you.

Image Analyst 2025-1-25

@sajjad, did you see my comment above: https://www.mathworks.com/matlabcentral/answers/36171-calculating-uaci-npcr-for-2-images#comment_991829

That's about all I can offer, is that File Exchange entry. I don't have any code for it myself.

请先登录，再进行评论。

Calculating UACI, NPCR for 2 images

5 个评论
显示 3更早的评论隐藏 3更早的评论

采纳的回答

2 个评论
显示无隐藏无

更多回答（1 个）

14 个评论
显示 12更早的评论隐藏 12更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Calculating UACI, NPCR for 2 images

5 个评论 显示 3更早的评论隐藏 3更早的评论

采纳的回答

2 个评论 显示 无隐藏 无

更多回答（1 个）

14 个评论 显示 12更早的评论隐藏 12更早的评论

另请参阅

类别

标签

Community Treasure Hunt

5 个评论
显示 3更早的评论隐藏 3更早的评论

2 个评论
显示无隐藏无

14 个评论
显示 12更早的评论隐藏 12更早的评论