fractal-mosaics/reference/upstream-matlab/fractal_mosaic.m
copilot 5bfeae238d Add Python fractal-mosaic CLI and unpack MATLAB reference
Fork of s-ben/Fractal-Mosaics for Tower mosaik inbox: GPL license,
upstream .m sources under reference/, RST placement with CBIR +
log-polar hints, edge alpha, profiles wenig/mittel/dicht, no recolor.

Co-authored-by: Cursor <cursoragent@cursor.com>
2026-09-24 17:03:12 +02:00

1127 lines
50 KiB
Matlab

% Copyright (c) 2013 Seth Benton (seth.benton@gmail.com)
%
% This program (Fractal Mosaics) is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program. If not, see <http://www.gnu.org/licenses/>.
% The latest source code and documentation for Fractal Mosaics can be
% found at: https://github.com/s-ben/Fractal-Mosaics
clear all
close all
load DC_vals_graffiti
tic
% ---------- User Parameters ----------------
% Use these parameters to control the behavior of the Fractal Mosaic
% algorithm
color = 0; % if color = '0', black and white images are used for matching. color = '1' uses color images
threshold_layer_on = 0; % '1' if a "transparency layer" (thresh_layer) is used to control the acceptable error (rms_bias_v) on a pixel-by-pixel basis.
% if threshold_layer_on = '0', the same error
target_path = ['rose_black_bg__corner_2048.jpg']; % path to target image
threshold_path = ['eye_pic_clean_thresh_layer.jpg']; % path to (optional) transparency image
% path = ['../library_images/san_francisco_256_color/']; % path to library images
path = ['../library_images/graffiti_256_black_white/']; % path to library images
bias_floor = 24; % bias floor
rms_bias_range = 6; % range of scaling according to transparency layer (0-1 greyscale image)
max_repeats_v = [1 2 4 4 20 50 100 100]; % maximum number of repeated images in the mosaic (level dependent)
keeper_thresh_v = [10 10 10 11 14 15 15 15]; % for each of the eight "levels" [1 2 3 4 5 6 7 8] (scales of analysis), set an rms threshold.
num_alt_images = 3; % number of alternate images to store per pixel
perc_max_overlap_v = [.20 .20 .20 .20 .30 .35 .40 .50]; % max overlap of candidate images on placed image
num_lib_test_v = [1 1 1 1 1 1 1 1]; % number of candidate matches to stop at at each level (i.e., if '1' for a given level, the algorithm won't continue to look
% "underneath" a placed pixel for more candidate images. You probably
% won't touch this parameter
num_lib_test_pix_v = [3000 2500 2000 1000 750 1000 500 500]; % number of library images to test for each pixel, per level
num_lib_test_pix_v = num_lib_test_pix_v/2; % REMOVE
% ------- resolution levels ------------
res_v = [8 16 32 64 128 256];
level_v = [5 5 5 5 5 5 4 3]; % for each level, the res_v() to go up to
% Create matrix for calculation of color histograms
range = 0.0:0.1:1.0;
rangeNew = 0.0:0.05:1.0;
[x1,y1,z1] = meshgrid(range);
[x2,y2,z2] = meshgrid(rangeNew);
% --------- rendering parameters -------------------
rendering_level = 5; % indexes res_v (e.g., 1 -> 8x8, 2 -> 16x16, ....)
scaling_constant_v = [7 8 9 10 11 12 13 14]; % vector of values used to scale candidate pixels to their respective analysis levels (j)
lowest_res_image = res_v(rendering_level);%128; % used below to calculate size of "analysis tiles" for each Level
% --------- calculate size of analysis tile for each level ---------------------
for i=1:length(level_v)
lib_im_v(i) = lowest_res_image*2^(length(level_v)-2-i);
end
% Define vector for storing rms error values
rms_err_v = zeros(length(max(num_lib_test_pix_v)),1);
% Calculate rms thresholds for each Level (i) and layer (j)
rms_thresh_v = [(bias_floor) (bias_floor) (bias_floor) (bias_floor) (bias_floor+1) ; % i=1
(bias_floor) (bias_floor) (bias_floor) (bias_floor) (bias_floor+1); % i=2
(bias_floor) (bias_floor) (bias_floor) (bias_floor) (bias_floor+1); % i=3
(bias_floor) (bias_floor) (bias_floor) (bias_floor) (bias_floor+3); % i=4
(bias_floor+2) (bias_floor+2) (bias_floor+2) (bias_floor+2) (bias_floor+4); % i=5
(bias_floor+2) (bias_floor+2) (bias_floor+2) (bias_floor+2) (bias_floor+6); % i=6
(bias_floor+4) (bias_floor+4) (bias_floor+4) (bias_floor+7) (bias_floor+7); % i=7
(bias_floor+7) (bias_floor+7) (bias_floor+10) (bias_floor+10) (bias_floor+10)]; % i=8
% number of candidate pixels (via level 'i' and resolution 'j') that are passed through to the next level
pass_through_v_abs = [1 4 4 1 1 1 ; 1 4 4 1 1 1 ; 1 4 4 1 1 1; 1 4 4 1 1 1; 1 4 1 1 1 1; 1 1 1 1 1 1; 1 1 1 1 1 1; 1 1 1 1 1 1];
% ------------- Load target image ------------------
if color == 0
I1 = imread(target_path);
I1_size = size(I1);
% If target an RGB image, converty to greyscale
if (I1_size(3) == 3)
I1 = rgb2gray(I1);
end
if threshold_layer_on == 1
thresh_layer = imread(threshold_path);
thresh_layer = rgb2gray(thresh_layer);
end
else
% I1 = imread('old_car_blue_2048_clean.jpg');
% thresh_layer = imread('old_car_blue_2048_clean_transparency.jpg');
I1 = imread(target_path);
I1_size = size(I1);
if threshold_layer_on == 1
thresh_layer = imread(threshold_path);
end
end
Im_target = I1;
clear I1;
% ----------------------------- get jpeg file paths --------------------
'getting image paths...'
directories = dir(path); % get directory names
num_files = 0;
file_count = 0;
dir_temp = directories(3).name; % get first directory name
files = dir([path,dir_temp,'/*.jpg']); % get file names for first directory
file_count = length(files);
% add directory path to jpeg files name
for j=1:length(files)
files(j).name = [path,dir_temp,'/',files(j).name];
end
for i=2:(length(directories)-2) % offset by 2 because first two entries are '.' and '..'
dir_temp = directories(i+2).name; % get next directory name
files_temp = dir([path,dir_temp,'/*.jpg']); % get file names for given directory
% add directory path to jpeg files name
for j=1:length(files_temp)
files_temp(j).name = [path,dir_temp,'/',files_temp(j).name];
end
% concatenate file paths from next directory to existing
files((file_count+1):(file_count+length(files_temp))) = files_temp;
file_count = file_count+length(files_temp);
end
% --- create varying resolution target images for matching ------
if color == 0
size_target = size(Im_target);
else
size_target = size(Im_target(:,:,1));
end
partitions = I1_size(1)/1024;
res_scale = (res_v(1)*partitions)/size_target(1); % height of lowest res target image for matching / orig target size
Im_target_lr = imresize(Im_target,res_scale); % lowest resolution target image
size_target_low_res = size(Im_target_lr);
for i=1:10 % for each resolution needed (hard-coded for now, need to generalize)
% ------- resize target to res -------
scale_factor = 2^(i-1);
Im_target_tmp(i).lr = imresize(Im_target,[size_target_low_res(1)*scale_factor size_target_low_res(2)*scale_factor]);
end
% ---------- create final 'Render_target' variable ------
Render_target = Im_target_tmp(10).lr;
size_render_target = size(Render_target);
% Render_target_valid is a logical bitmap specifying whether an images has been "placed" on the target image.
% declar all '0s' logical image
Render_target_valid = im2bw(Render_target,1);
% ---------------------------------------------------------------
if color ==0
% -------- load histograms of library images ------------
load('hist_matrix_graffiti_128');
% load('hist_matrix_SF_128');
% ------------- Load pdfs for hist_dist, gfd, std_diff ------------------
load('pdfs');
% -------- load GFDs of library images ------------
load('gfd_matrix_graffiti_128');
% load('gfd_matrix_SF_128');
load('std_matrix_graffiti_128');
% load('std_matrix_SF_128');
else
load san_francisco_color_hists;
%load graffiti_color_hists;
end
% -------- resize threshold layer (assuming i=8 is highest) ----------------
Im_size_lr = size(Im_target_tmp(8).lr);
if color == 0
if threshold_layer_on == 1
thresh_layer = imresize(thresh_layer, Im_size_lr);
end
else
if threshold_layer_on == 1
thresh_layer = imresize(thresh_layer, [Im_size_lr(1) Im_size_lr(2)]);
end
end
% ---------- calculate distance windows -----------
for i=1:length(level_v)
perc_sloped = .05;
perc_sloped_fixed = .01; % perc sloped for fixed window (one that's used to cut off clipping)
fd = zeros(lib_im_v(i));
fd_fixed = zeros(lib_im_v(i));
y_d = ones(1,lib_im_v(i)/2);
y_d_fixed = ones(1,lib_im_v(i)/2);
x_ind = linspace(0,lib_im_v(i)/2,lib_im_v(i)/2);
for u=1:(perc_sloped*(lib_im_v(i)/2))
y_d(u) = (1/(perc_sloped*(lib_im_v(i)/2)))*x_ind(u);
end
for u=1:(perc_sloped_fixed*(lib_im_v(i)/2))
y_d_fixed(u) = (1/(perc_sloped_fixed*(lib_im_v(i)/2)))*x_ind(u);
end
fd(1:lib_im_v(i)/2) = y_d;
fd(lib_im_v(i)/2+1:lib_im_v(i)) = fliplr(y_d);
fd_2 = fd*fd';
window(i).win = fd_2;
fd_fixed(1:lib_im_v(i)/2) = y_d_fixed;
fd_fixed(lib_im_v(i)/2+1:lib_im_v(i)) = fliplr(y_d_fixed);
fd_2_fixed = fd_fixed*fd_fixed';
window(i).win_fixed = fd_2_fixed;
end
% ----- allocate array to store candidate pixels -------
cand_pixels = zeros(2,4^(length(res_v)-1));
cand_pixels_next = zeros(2,4^(length(res_v)-1));
cand_pixels_rms = zeros(2,4^(length(res_v)-1));
'rendering mosaic'
% set counters to zero
num_cand_pixels_final = 0;
num_cand_pixels_alt = 0;
cand_pixels_final_abs = zeros(10000,3);
dist = zeros(1,30000);
% calculate histogram bins
bins = 16;
hist_x = [(255/bins)/2:255/bins:(255-(255/bins)/2)];
% declare final pdf (sum of pdfs) array
match_pdf = zeros(length(files),1);
% declare sq_center for each level (i.e., declare pixel boundaries)
% the pixel boundaries (for now) manually squeeze in to avoid too many
% lower level (time consuming) matches on the borders
sq_center_i_v = [5 5 5 7 14 28 56 112];
% vector of scalings for rad_buff (establishes radius for palced matches)
buff_v = [3 3 3.5 3.8 3.9 4 4.2 4.5];
for i=2:length(level_v) % for each scale (level)
% 'scale'
% i
% ------- read lowest resolution target image for given scale -------
Im_target_match = Im_target_tmp(i).lr;
if color ==0
size_Im_target_match = size(Im_target_match);
else
size_Im_target_match = size(Im_target_match(:,:,1));
end
% ------- read in target image w/ 128 resolution for current scale ------
Im_target_match_hist = Im_target_tmp((i-1)+level_v(i)).lr;
% ----- for each pixel in valid area ---------
for y = sq_center_i_v(i):(size_Im_target_match(1)-sq_center_i_v(i)+1) % for each "translation" pixel in (lowest res) valid area (y_sample(v)*(2^(i-1))):(y_sample(v)*(2^(i-1)))
for x = sq_center_i_v(i):(size_Im_target_match(2)-sq_center_i_v(i)+1) %sq_center_i(2):(size_Im_target_match(2)-sq_center_i(2)+1)
img_tested = num_lib_test_v(i); % number of library images at this level to evaluate
num_img_in_rad = 0; % number of images within the radius of already selected candidate images
quota_reached = 0; % binary variable to flag if quota has been reached
cand_image_found = 0;
num_cand_image_found = 0; % number of candidate images found
% ---------- scale cand pixel position according to scale and res ---------------------
x_tmper = round(x*2^(scaling_constant_v(rendering_level)-(i+1)));
y_tmper = round(y*2^(scaling_constant_v(rendering_level)-(i+1)));
% ------------ find rms bias for given pixel ------------------------
x_thresh = x*2^(8-i);
y_thresh = y*2^(8-i);
% -------------------------------------------------------------
% If threshold layer turned off, don't vary rms error threshold
% pixel-by-pixel
if (threshold_layer_on == 0)
rms_bias_v = 0;
else
rms_bias_v = thresh_layer(y_thresh,x_thresh);
rms_bias_v = rms_bias_range*(1-(double(rms_bias_v)/255));
end
% ------------- calculate distance of cand pixel to placed images --------------------------------
for w=1:num_cand_pixels_final
% If a library image has already been placed which covers
% the candidate pixel, exit loop
if (Render_target_valid(y_tmper,x_tmper) == 1)
quota_reached = 1;
break
end
% calculate the minimum distance from the center of a place
% image that a candidate image needs to be to be evaluated
min_dist = cand_pixels_final_abs(w,3) - perc_max_overlap_v(i)*lib_im_v(i) + lib_im_v(i)/2;
% calculate distance from candidate pixel location to
% placed image location
dist(w) = sqrt((cand_pixels_final_abs(w,1) - x_tmper)^2 + (cand_pixels_final_abs(w,2) - y_tmper)^2);
% if candidate pixel is too close to a placed image, exit
% loop
if (dist(w) < min_dist)
num_img_in_rad = num_img_in_rad+1;
if (num_img_in_rad >= num_lib_test_v(i))
quota_reached = 1;
break
end
end
end
% --------------- if quota reached for cand pixels in this area, skip loop (don't evaluate this cand pixel) -------------------------
if (quota_reached == 1)
% 'quota reached, exiting loop'
continue
% z = 1000000;
else
img_tested = num_lib_test_v(i)-num_img_in_rad;
z = 1;
end
% ------- calculate 128x128 target_tile for current pixel --------
x_tmp = x;
y_tmp = y;
for z=2:(level_v(i))
% ----update pixel location for higher res ---------
x_tmp = 2*(x_tmp-0.5);
y_tmp = 2*(y_tmp-0.5);
end
xmin = x_tmp-(floor((res_v(level_v(i)))/2)+1) +1;
ymin = y_tmp-(floor((res_v(level_v(i)))/2)+1) +1;
width = res_v(level_v(i))-1;
height = res_v(level_v(i))-1;
rect = [xmin ymin width height];
Im_target_tile_hist = imcrop(Im_target_match_hist,rect);
if (res_v(level_v(i))< res_v(rendering_level))
Im_target_tile_hist = imresize(Im_target_tile_hist,[res_v(rendering_level) res_v(rendering_level)]);
end
% -------------------------------------------------------------
% ------------- calculate hist, std, and GFD of 'target_tile' -----------------------
% Parameters for calculating General Fourier Descriptor (GFD)
radial_freq = 6;
angular_freq = 6;
size_v = [128 128];
if color==0
% -------- calculate DC component --------------
DC_component = sum(sum(Im_target_tile_hist))/(128*128);
% ----- calculate histogram ----------
Im_lib_array = reshape(Im_target_tile_hist,128*128,1);
hist_target_tile = hist(double(Im_lib_array),hist_x);
if (max(max(Im_target_tile_hist)) == 0)
Im_target_tile_hist = Im_target_tile_hist + 0.01;
end
% ------ calculate GFD -------------
Im_lib_polar = polar_conversion(Im_target_tile_hist);
Im_lib_polar_FFT = abs(fftshift(fft2(Im_lib_polar)));
DC = Im_lib_polar_FFT(size_v(1)/2+1,size_v(2)/2+1);
Im_lib_polar_FFT = Im_lib_polar_FFT./DC;
center_gfd = [size_v(1)/2+1 size_v(2)/2+1];
gfd = Im_lib_polar_FFT(center_gfd(1):(center_gfd(1)+angular_freq-1),center_gfd(2):(center_gfd(2)+radial_freq-1));
% ------ calculate std -------------
std_tile = std2(Im_target_tile_hist);
else
% -------- Get color image histogram ----------------
RGB = Im_target_tile_hist;
RGBt = RGB;
RGB = rgb2hsv(RGB);
% get image size:
[M,N,ttt] = size(RGB);
range = 0.0:0.1:1.0;
Hist = zeros(length(range),length(range),length(range));
for ii=1:M
for jj=1:N
nn1 = round(RGB(ii,jj,1) * 10)+1;
nn2 = round(RGB(ii,jj,2) * 10)+1;
nn3 = round(RGB(ii,jj,3) * 10)+1;
Hist(nn1, nn2, nn3) = Hist(nn1, nn2, nn3) + 1;
end
end
Hist = Hist / (M*N);
% Calculate DC component
img_tmp = rgb2gray(Im_target_tile_hist);
DC_component = sum(sum(img_tmp))/(256*256);
% ----------------------------------------------------------
end
% ------------ calculate statistical difference (hist, gfd, std) ---------------
if color ==0
% ------------ calculate histogram distance ---------------
neighbor_distances = kNearestNeighbors_no_sort(hist_matrix,hist_target_tile,length(files));
% ------------ calculate GFD distance ---------------
neighbor_distances_gfd = kNearestNeighbors_no_sort(gfd_matrix,gfd(:)',length(files));
% ------------ calculate std distance ---------------
neighbor_distances_std = abs(std_tile - std_matrix);
% neighbor_distances_int = neighbor_distances;
else
Nfiles = length(Hists);
% decision thresholds:
t = 0.010;
t2 = 0.8;
Hist = interp3(x1,y1,z1,Hist,x2,y2,z2);
Similarity = zeros(Nfiles, 1);
for (w=1:Nfiles) % for each file in database:
% compute (normalized) eucledean distance for all hist bins:
HistT = Hists{w};
% HistT = interp3(x,y,z,Hists{i},x2,y2,z2);
DIFF = abs(Hist-HistT) ./ Hist;
% keep distance values for which the corresponding query image's values
% are larger than the predefined threshold:
DIFF = DIFF(Hist>t);
% keep error values which are smaller than 1:
DIFF2 = DIFF(DIFF<t2);
L2 = length(DIFF2);
% compute the similarity meaasure:
Similarity(w) = length(DIFF) * mean(DIFF2) / (L2^2);
end
% find the nResult "closest" images:
% [Sorted, ISorted] = sort(Similarity);
end
% ---------------- evaluate match probability from stats -------------------------------
if color==0
for z=1:length(hist_dist_bin)
hist_array(:,z) = abs(neighbor_distances - hist_dist_bin(z));
end
for z=1:length(gfd_bin)
gfd_array(:,z) = abs(neighbor_distances_gfd - gfd_bin(z));
end
for z=1:length(std_diff_bin)
std_array(:,z) = abs(neighbor_distances_std - std_diff_bin(z));
end
% Find which histogram bin the measured statistic (histogram,
% gfd, std) are nearest to
[hist_min hist_ind] = min(hist_array,[],2);
[gfd_min gfd_ind] = min(gfd_array,[],2);
[std_min std_ind] = min(std_array,[],2);
match_pdf = hist_dist_pdf(hist_ind) + gfd_pdf(gfd_ind) + std_diff_pdf(std_ind);
end
if color==0
[sortval nearest_neighbor_int_index] = sort(match_pdf,'descend');
else
[sortval nearest_neighbor_int_index] = sort(Similarity);
% [sortval nearest_neighbor_int_index] = sort(neighbor_distances);
end
max_rms_final = 1000; % reset max rms for final resolution (if pixel reaches it)
z=1;
% ----------- Calculate 8x8 registration (retrieval step) -------------------------
j=1;
% ------- read in target image w/ correct resolution ------
Im_target_match = Im_target_tmp((i-1)+j).lr;
% ------ calculate valid pixel area -------
sq_center = [floor(res_v(j)/2)+1 floor(res_v(j)/2)+1];
j=1; % only temporarily. gets reassigned in for loop below
xmin = x-sq_center(2) +1;
ymin = y-sq_center(1) +1;
width = res_v(j)-1;
height = res_v(j)-1;
rect = [xmin ymin width height];
% crop piece of target image we're looking to match to
Im_target_tile = imcrop(Im_target_match,rect);
% Take log-polar transform of Im_target_tile
if color==0
L1 = transformImage_log(Im_target_tile, res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
else
L1_r = transformImage_log(Im_target_tile(:,:,1), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
L1_g = transformImage_log(Im_target_tile(:,:,2), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
L1_b = transformImage_log(Im_target_tile(:,:,3), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
end
% ----------- set up scale finding function ----------
minScale = 1;
maxScale = res_v(j)/2-1;
rad_lin = linspace(log(minScale), log(maxScale),res_v(j))';
rho = exp(rad_lin);
[min_difference, array_position] = min(abs(rho - 2));
Im_target_tile_masked = zeros([res_v(j) res_v(j)]);
z=1;
% for each library image tested, calculated rms error
while(z<=num_lib_test_pix_v(i))
% calculate number of repeated images at the current level
repeats = 0;
if (nearest_neighbor_int_index(z) == 9522) % if z points to a "bad image" (e.g. found to be boring), skip
z = z+1;
else % check for repeat images
for p=1:num_cand_pixels_final
if (z == cand_pixels_final(num_cand_pixels_final).lib_im)
repeats = repeats +1;
end
if (repeats >= max_repeats_v(i))
z=z+1;
end
end
end
image_path = [files(nearest_neighbor_int_index(z)).name];
for j=1:length(image_path)
if (image_path(j) == '\')
image_path(j) = '/';
end
end
Im_lib = imread(image_path); % read in library image
perc_DC = 1; % percentage DC adjusted
DC_adjust = (DC_v(nearest_neighbor_int_index(z))-DC_component)*perc_DC;
Im_lib = Im_lib-DC_adjust;
for j=1:level_v(i) % for each resolution (up to the highest acceptable for each scale)
% ------- read in target image w/ correct resolution ------
Im_target_match = Im_target_tmp((i-1)+j).lr;
% -------resize library image-------------
Im_lib_tmp = imresize(Im_lib, [res_v(j) res_v(j)]);
% Im_lib_tmp = Im_lib;
% ------ calculate valid pixel area -------
sq_center = [floor(res_v(j)/2)+1 floor(res_v(j)/2)+1];
% ------ get size of current res target image -------
if color == 0
size_target = size(Im_target_match);
else
size_target = size(Im_target_match(:,:,1));
end
% calculate the number of pixels to evaluate at this
% resolution
if j==1
num_cand_pixels = 1;
else
[rms_sort rms_index] = sort(cand_pixels_rms(1:num_cand_pixels_new));
num_cand_pixels = min(num_cand_pixels_new, pass_through_v_abs(i,j));
end
if num_cand_pixels == 0
break
end
% ---------- assign candidates from last pass for anaylsis ------------
if j>1
for b=1:num_cand_pixels
% cand_pixels(1,b) = cand_pixels_next(1,b);
% cand_pixels(2,b) = cand_pixels_next(2,b);
cand_pixels(1,b) = cand_pixels_next(1,rms_index(b));
cand_pixels(2,b) = cand_pixels_next(2,rms_index(b));
end
end
% % ------------- calculate log-polar of target tile ----------------------------
% ---- log polar transform library image ------
if color==0
L2 = transformImage_log(Im_lib_tmp, res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
else
L2_r = transformImage_log(Im_lib_tmp(:,:,1), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
L2_g = transformImage_log(Im_lib_tmp(:,:,2), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
L2_b = transformImage_log(Im_lib_tmp(:,:,3), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
end
num_cand_pixels_new = 0;
for b=1:num_cand_pixels
% ---------- if desired number of candidate images have been reached, break from search -----------------------
% ----- update pixel position at higher resolutions (if necessary) -----
if j==1
x_c = x;
y_c = y;
else
%--- read in candidate pixel location---
y_c = cand_pixels(1,b);
x_c = cand_pixels(2,b);
% ----update pixel location for higher res ---------
x_c = 2*(x_c-0.5);
y_c = 2*(y_c-0.5);
end
% ------ determine indexing parameters for higher resolution (if necessary) ------
% res_scale_up = res_v(j)/res_v(1);
if j==1
x_start = x_c;
x_end = x_c;
y_start = y_c;
y_end = y_c;
pix_dist = 1;
else
pix_dist = 2;
x_start = x_c;
x_end = x_start + pix_dist-1;
y_start = y_c;
y_end = y_start + pix_dist-1;
end
% ----- for each sub-pixel (if not lowest res) --------------
for p = y_start:y_end
for q = x_start:x_end
if color==0
Im_target_tile_masked = zeros([res_v(j) res_v(j)]);
else
Im_target_tile_masked = zeros([res_v(j) res_v(j) 3]);
end
valid_pix = 0;
mask = zeros([res_v(j) res_v(j)]);
xmin = q-sq_center(2) +1;
ymin = p-sq_center(1) +1;
width = res_v(j)-1;
height = res_v(j)-1;
rect = [xmin ymin width height];
Im_target_tile = imcrop(Im_target_match,rect);
if color==0
L1 = transformImage_log(Im_target_tile, res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
else
L1_r = transformImage_log(Im_target_tile(:,:,1), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
L1_g = transformImage_log(Im_target_tile(:,:,2), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
L1_b = transformImage_log(Im_target_tile(:,:,3), res_v(j), res_v(j), res_v(j), res_v(j), 'bilinear', [res_v(j) res_v(j)] / 2, 'valid');
end
if color==0
cc=ifft2(fft2(L1).*conj(fft2(L2)));
cc=flipud(cc);
cc = circshift(cc,[(1) (0)]);
else
cc=ifft2(fft2(L1_r).*conj(fft2(L2_r)))+ifft2(fft2(L1_g).*conj(fft2(L2_g)))+ifft2(fft2(L1_b).*conj(fft2(L2_b)));
cc=flipud(cc);
cc = circshift(cc,[(1) (0)]);
end
% ----------- set up scale finding function ----------
minScale = 1;
maxScale = res_v(j)/2-1;
rad_lin = linspace(log(minScale), log(maxScale),res_v(j))';
rho = exp(rad_lin);
[min_difference, array_position] = min(abs(rho - 2));
% -------- find max of cross correlation ---------------------
[y_max, x_ind] = max(cc(1:array_position,:)); % returns max values per row along with x-index value
[x_max, y_ind] = max(y_max);
y_max_index = x_ind(y_ind);
x_max_index = y_ind;
% --------- calculate rotation degree -------------------------
Theta = 360*((x_max_index -1)/res_v(j));
% --------- calculate scale value ----------------------------
scale_est = rho(y_max_index);
I_sc_tmp = imresize(Im_lib_tmp, 1/scale_est);
I_sc_tmp_size = size(I_sc_tmp);
tmp_height = I_sc_tmp_size(1);
tmp_width = I_sc_tmp_size(2);
mask_tmp = ones([I_sc_tmp_size(1) I_sc_tmp_size(2)]);
if color==0
I_sc = zeros([res_v(j) res_v(j)]);
I_comb = zeros([res_v(j) res_v(j)]);
else
I_sc = zeros([res_v(j) res_v(j) 3]);
I_comb = zeros([res_v(j) res_v(j) 3]);
end
if color==0
I_sc(round(res_v(j)/2 - tmp_height/2)+1:round(res_v(j)/2 - tmp_height/2)+tmp_height,round(res_v(j)/2 - tmp_width/2)+1:round(res_v(j)/2 - tmp_width/2)+tmp_width) = I_sc_tmp;
else
I_sc(round(res_v(j)/2 - tmp_height/2)+1:round(res_v(j)/2 - tmp_height/2)+tmp_height,round(res_v(j)/2 - tmp_width/2)+1:round(res_v(j)/2 - tmp_width/2)+tmp_width,:) = I_sc_tmp;
end
mask(round(res_v(j)/2 - tmp_height/2)+1:round(res_v(j)/2 - tmp_height/2)+tmp_height,round(res_v(j)/2 - tmp_width/2)+1:round(res_v(j)/2 - tmp_width/2)+tmp_width) = mask_tmp;
% -------Rotate image back by theta estimate --------
R1 = imrotate(I_sc, -Theta, 'bilinear', 'crop');
mask = imrotate(mask, -Theta, 'bilinear','crop');
alpha = 0.5;
I_comb = (alpha)*double(Im_target_tile) + (1-alpha)*R1;
if color==0
Im_target_tile_masked = double(Im_target_tile).*mask;
else
Im_target_tile_masked(:,:,1) = double(Im_target_tile(:,:,1)).*mask;
Im_target_tile_masked(:,:,2) = double(Im_target_tile(:,:,2)).*mask;
Im_target_tile_masked(:,:,3) = double(Im_target_tile(:,:,3)).*mask;
end
rms_err = sqrt( (sum(sum( (abs(R1-Im_target_tile_masked)).^2 )))/(sum(sum(mask))));
if color==1
rms_err = sum(rms_err)/3;
end
rms_err_v(z,j) = rms_err;
if rms_err < (rms_thresh_v(i,j)-rms_bias_v)
num_cand_pixels_new = num_cand_pixels_new+1;
% -------store candidate pixels for next pass ---------
cand_pixels_next(1,num_cand_pixels_new) = p;
cand_pixels_next(2,num_cand_pixels_new) = q;
cand_pixels_rms(num_cand_pixels_new) = rms_err;
% ---------- If candidate pixel is in highest res for given scale, save ------------
if (j == level_v(i))
% ---------- scale cand pixel position according to scale and res ---------------------
x_tmper = round(q*2^(scaling_constant_v(rendering_level)-(i+j)));
y_tmper = round(p*2^(scaling_constant_v(rendering_level)-(i+j)));
radius = ((lib_im_v(i)/2)/scale_est);
if rms_err < max_rms_final
% ----------- if >1 candidates, store alternatives--------------
if ((num_cand_image_found > 1)&& (num_cand_image_found <= num_alt_images))
num_cand_pixels_alt = num_cand_pixels_alt+1;
cand_pixels_final_alt(num_cand_pixels_alt).DC = cand_pixels_final(num_cand_pixels_final+1).DC;
cand_pixels_final_alt(num_cand_pixels_alt).rms = cand_pixels_final(num_cand_pixels_final+1).rms;
cand_pixels_final_alt(num_cand_pixels_alt).x_cand = cand_pixels_final(num_cand_pixels_final+1).x_cand;
cand_pixels_final_alt(num_cand_pixels_alt).y_cand = cand_pixels_final(num_cand_pixels_final+1).y_cand;
cand_pixels_final_alt(num_cand_pixels_alt).res = cand_pixels_final(num_cand_pixels_final+1).res;
cand_pixels_final_alt(num_cand_pixels_alt).scale = cand_pixels_final(num_cand_pixels_final+1).scale;
cand_pixels_final_alt(num_cand_pixels_alt).rot = cand_pixels_final(num_cand_pixels_final+1).rot;
cand_pixels_final_alt(num_cand_pixels_alt).lib_im = cand_pixels_final(num_cand_pixels_final+1).lib_im;
cand_pixels_final_alt(num_cand_pixels_alt).level = cand_pixels_final(num_cand_pixels_final+1).level;
cand_pixels_final_alt(num_cand_pixels_alt).valid_pix = cand_pixels_final(num_cand_pixels_final+1).valid_pix;
cand_pixels_final_alt(num_cand_pixels_alt).file_path = cand_pixels_final(num_cand_pixels_final+1).file_path;
end
cand_pixels_final(num_cand_pixels_final+1).DC = DC_adjust;
cand_pixels_final(num_cand_pixels_final+1).rms = rms_err;
cand_pixels_final(num_cand_pixels_final+1).x_cand = q;
cand_pixels_final(num_cand_pixels_final+1).y_cand = p;
cand_pixels_final(num_cand_pixels_final+1).res = j;
cand_pixels_final(num_cand_pixels_final+1).scale = scale_est;
cand_pixels_final(num_cand_pixels_final+1).rot = Theta;
cand_pixels_final(num_cand_pixels_final+1).lib_im = nearest_neighbor_int_index(z);
cand_pixels_final(num_cand_pixels_final+1).level = i;
cand_pixels_final(num_cand_pixels_final+1).valid_pix = valid_pix;
cand_path = files(nearest_neighbor_int_index(z)).name;
cand_pixels_final(num_cand_pixels_final+1).file_path = cand_path(length(path):end);
max_rms_final = rms_err;
cand_image_found = 1;
num_cand_image_found = num_cand_image_found+1; % count number of images below rms threshold for given pixel
cand_pixels_final_abs(num_cand_pixels_final+1,1) = x_tmper;
cand_pixels_final_abs(num_cand_pixels_final+1,2) = y_tmper;
cand_pixels_final_abs(num_cand_pixels_final+1,3) = radius;
end
end
end
% figure(1),subplot(2,2,1),imshow(Im_lib_tmp), title('Im_lib')
% subplot(2,2,2),imshow(Im_target_tile), title('Target Tile')
% subplot(2,2,3),imshow(uint8(R1)), title('transformed Image')
% subplot(2,2,4),imshow(uint8(I_comb)), title('Combined Images')
% pause
%
% rms_err
end
end
end
end
z = z+1;
end
% ------------ if highest res math found, increment final pixel count, Render match -------------
if (cand_image_found ==1)
num_cand_pixels_final = num_cand_pixels_final+1;
% ------------ Render match ----------------
DC_adjust = cand_pixels_final(num_cand_pixels_final).DC;
i_r = cand_pixels_final(num_cand_pixels_final).level;
j_r = cand_pixels_final(num_cand_pixels_final).res;
x_r = cand_pixels_final(num_cand_pixels_final).x_cand;
y_r = cand_pixels_final(num_cand_pixels_final).y_cand;
scale_est_r = cand_pixels_final(num_cand_pixels_final).scale;
Theta_r = cand_pixels_final(num_cand_pixels_final).rot;
lib_im_num_r = cand_pixels_final(num_cand_pixels_final).lib_im;
% ---------- scale cand pixel position according to scale and res ---------------------
x_r = round(x_r*2^(scaling_constant_v(rendering_level)-(i_r+j_r)));
y_r = round(y_r*2^(scaling_constant_v(rendering_level)-(i_r+j_r)));
% ---------- read in candidate library image ------------------
Im_lib = imread([files(lib_im_num_r).name]); % read in library image
perc_DC = 1; % percentage DC adjusted
% DC_adjust = (DC_v(nearest_neighbor_int_index(z))-DC_component)*perc_DC;
Im_lib = Im_lib-DC_adjust;
% -------resize library image for level -------------
Im_lib_tmp = imresize(Im_lib, [lib_im_v(i_r) lib_im_v(i_r)]);
% ------- crop analysis window -------------------------
xmin = x_r - lib_im_v(i_r)/2;
ymin = y_r - lib_im_v(i_r)/2;
width = lib_im_v(i_r)-1;
height = lib_im_v(i_r)-1;
rect = [xmin ymin width height];
Im_target_tile = imcrop(Render_target,rect);
Im_target_tile_masked = zeros([lib_im_v(i_r) lib_im_v(i_r)]);
mask = zeros([lib_im_v(i_r) lib_im_v(i_r)]);
% ---------- matrix for backing up Render target pixels -----------------
Render_target_backup = zeros([lib_im_v(i_r) lib_im_v(i_r)]);
% --------- scale image ----------------
I_sc_tmp = imresize(Im_lib_tmp, 1/scale_est_r);
I_sc_tmp_size = size(I_sc_tmp);
tmp_height = I_sc_tmp_size(1);
tmp_width = I_sc_tmp_size(2);
radius = ((lib_im_v(i_r)/2)/scale_est_r)*(1-perc_sloped);
fd_2 = window(i).win;
if color == 0
fd_2 = imresize(fd_2,I_sc_tmp_size);
else
fd_2 = imresize(fd_2,[I_sc_tmp_size(1) I_sc_tmp_size(2)]);
end
if color == 0
mask_tmp = ones(I_sc_tmp_size);
else
mask_tmp = ones([I_sc_tmp_size(1) I_sc_tmp_size(2)]);
end
if color == 0
I_sc = zeros([lib_im_v(i_r) lib_im_v(i_r)]);
else
I_sc = zeros([lib_im_v(i_r) lib_im_v(i_r) 3]);
end
I_comb = zeros([lib_im_v(i_r) lib_im_v(i_r)]);
if color == 0
I_sc(round(lib_im_v(i_r)/2 - tmp_height/2)+1:round(lib_im_v(i_r)/2 - tmp_height/2)+tmp_height,round(lib_im_v(i_r)/2 - tmp_width/2)+1:round(lib_im_v(i_r)/2 - tmp_width/2)+tmp_width) = I_sc_tmp;
else
I_sc(round(lib_im_v(i_r)/2 - tmp_height/2)+1:round(lib_im_v(i_r)/2 - tmp_height/2)+tmp_height,round(lib_im_v(i_r)/2 - tmp_width/2)+1:round(lib_im_v(i_r)/2 - tmp_width/2)+tmp_width, :) = I_sc_tmp;
end
mask(round(lib_im_v(i_r)/2 - tmp_height/2)+1:round(lib_im_v(i_r)/2 - tmp_height/2)+tmp_height,round(lib_im_v(i_r)/2 - tmp_width/2)+1:round(lib_im_v(i_r)/2 - tmp_width/2)+tmp_width) = fd_2;
% -------Rotate image back by theta estimate --------
R1 = imrotate(I_sc, -Theta_r, 'bilinear', 'crop');
mask = imrotate(mask, -Theta_r, 'bilinear', 'crop');
fd_2_tmp = window(i).win_fixed;
mask = mask.*fd_2_tmp;
alpha = 0.5;
I_comb = (alpha)*double(Im_target_tile) + (1-alpha)*R1;
% ----------------- Render match ---------------------
for u= 1:lib_im_v(i_r)
for v = 1:lib_im_v(i_r)
Render_target_backup(u,v) = Render_target(ymin+u-1,xmin+v-1);
if (mask(u,v) > .01)
if color == 0
Render_target(ymin+u-1,xmin+v-1) = uint8(R1(u,v)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1));
else
Render_target(ymin+u-1,xmin+v-1,1) = uint8(R1(u,v,1)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1,1));
Render_target(ymin+u-1,xmin+v-1,2) = uint8(R1(u,v,2)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1,2));
Render_target(ymin+u-1,xmin+v-1,3) = uint8(R1(u,v,3)*mask(u,v) + (1-mask(u,v))*Render_target(ymin+u-1,xmin+v-1,3));
end
% update bitmap of valid placed pixels
Render_target_valid(ymin+u-1,xmin+v-1) = 1;
end
end
end
end
end
end
level_txt = 1;%num2str(i);
if (num_cand_pixels_final>0)
save(['Fractal_Mosaic_' level_txt])
end
end
% end
imwrite(Render_target,'fractal_mosaic','jpg','quality',100);
time = toc