diff --git a/imresize.py b/imresize.py new file mode 100644 index 0000000..bb88490 --- /dev/null +++ b/imresize.py @@ -0,0 +1,146 @@ +from __future__ import print_function +import numpy as np +from math import ceil, floor + +def deriveSizeFromScale(img_shape, scale): + output_shape = [] + for k in range(2): + output_shape.append(int(ceil(scale[k] * img_shape[k]))) + return output_shape + +def deriveScaleFromSize(img_shape_in, img_shape_out): + scale = [] + for k in range(2): + scale.append(1.0 * img_shape_out[k] / img_shape_in[k]) + return scale + +def triangle(x): + x = np.array(x).astype(np.float64) + lessthanzero = np.logical_and((x>=-1),x<0) + greaterthanzero = np.logical_and((x<=1),x>=0) + f = np.multiply((x+1),lessthanzero) + np.multiply((1-x),greaterthanzero) + return f + +def cubic(x): + x = np.array(x).astype(np.float64) + absx = np.absolute(x) + absx2 = np.multiply(absx, absx) + absx3 = np.multiply(absx2, absx) + f = np.multiply(1.5*absx3 - 2.5*absx2 + 1, absx <= 1) + np.multiply(-0.5*absx3 + 2.5*absx2 - 4*absx + 2, (1 < absx) & (absx <= 2)) + return f + +def contributions(in_length, out_length, scale, kernel, k_width): + if scale < 1: + h = lambda x: scale * kernel(scale * x) + kernel_width = 1.0 * k_width / scale + else: + h = kernel + kernel_width = k_width + x = np.arange(1, out_length+1).astype(np.float64) + u = x / scale + 0.5 * (1 - 1 / scale) + left = np.floor(u - kernel_width / 2) + P = int(ceil(kernel_width)) + 2 + ind = np.expand_dims(left, axis=1) + np.arange(P) - 1 # -1 because indexing from 0 + indices = ind.astype(np.int32) + weights = h(np.expand_dims(u, axis=1) - indices - 1) # -1 because indexing from 0 + weights = np.divide(weights, np.expand_dims(np.sum(weights, axis=1), axis=1)) + aux = np.concatenate((np.arange(in_length), np.arange(in_length - 1, -1, step=-1))).astype(np.int32) + indices = aux[np.mod(indices, aux.size)] + ind2store = np.nonzero(np.any(weights, axis=0)) + weights = weights[:, ind2store] + indices = indices[:, ind2store] + return weights, indices + +def imresizemex(inimg, weights, indices, dim): + in_shape = inimg.shape + w_shape = weights.shape + out_shape = list(in_shape) + out_shape[dim] = w_shape[0] + outimg = np.zeros(out_shape) + if dim == 0: + for i_img in range(in_shape[1]): + for i_w in range(w_shape[0]): + w = weights[i_w, :] + ind = indices[i_w, :] + im_slice = inimg[ind, i_img].astype(np.float64) + outimg[i_w, i_img] = np.sum(np.multiply(np.squeeze(im_slice, axis=0), w.T), axis=0) + elif dim == 1: + for i_img in range(in_shape[0]): + for i_w in range(w_shape[0]): + w = weights[i_w, :] + ind = indices[i_w, :] + im_slice = inimg[i_img, ind].astype(np.float64) + outimg[i_img, i_w] = np.sum(np.multiply(np.squeeze(im_slice, axis=0), w.T), axis=0) + if inimg.dtype == np.uint8: + outimg = np.clip(outimg, 0, 255) + return np.around(outimg).astype(np.uint8) + else: + return outimg + +def imresizevec(inimg, weights, indices, dim): + wshape = weights.shape + if dim == 0: + weights = weights.reshape((wshape[0], wshape[2], 1, 1)) + outimg = np.sum(weights*((inimg[indices].squeeze(axis=1)).astype(np.float64)), axis=1) + elif dim == 1: + weights = weights.reshape((1, wshape[0], wshape[2], 1)) + outimg = np.sum(weights*((inimg[:, indices].squeeze(axis=2)).astype(np.float64)), axis=2) + if inimg.dtype == np.uint8: + outimg = np.clip(outimg, 0, 255) + return np.around(outimg).astype(np.uint8) + else: + return outimg + +def resizeAlongDim(A, dim, weights, indices, mode="vec"): + if mode == "org": + out = imresizemex(A, weights, indices, dim) + else: + out = imresizevec(A, weights, indices, dim) + return out + +def imresize(I, scalar_scale=None, method='bicubic', output_shape=None, mode="vec"): + if method == 'bicubic': + kernel = cubic + elif method == 'bilinear': + kernel = triangle + else: + raise ValueError('unidentified kernel method supplied') + + kernel_width = 4.0 + # Fill scale and output_size + if scalar_scale is not None and output_shape is not None: + raise ValueError('either scalar_scale OR output_shape should be defined') + if scalar_scale is not None: + scalar_scale = float(scalar_scale) + scale = [scalar_scale, scalar_scale] + output_size = deriveSizeFromScale(I.shape, scale) + elif output_shape is not None: + scale = deriveScaleFromSize(I.shape, output_shape) + output_size = list(output_shape) + else: + raise ValueError('either scalar_scale OR output_shape should be defined') + scale_np = np.array(scale) + order = np.argsort(scale_np) + weights = [] + indices = [] + for k in range(2): + w, ind = contributions(I.shape[k], output_size[k], scale[k], kernel, kernel_width) + weights.append(w) + indices.append(ind) + B = np.copy(I) + flag2D = False + if B.ndim == 2: + B = np.expand_dims(B, axis=2) + flag2D = True + for k in range(2): + dim = order[k] + B = resizeAlongDim(B, dim, weights[dim], indices[dim], mode) + if flag2D: + B = np.squeeze(B, axis=2) + return B + +def convertDouble2Byte(I): + B = np.clip(I, 0.0, 1.0) + B = 255*B + return np.around(B).astype(np.uint8) + diff --git a/reweightedOptimizationSR.m b/reweightedOptimizationSR.m new file mode 100644 index 0000000..bede0b8 --- /dev/null +++ b/reweightedOptimizationSR.m @@ -0,0 +1,397 @@ +function [SR, model, report] = reweightedOptimizationSR(LRImages, model, optimParams, varargin) + + if nargin < 3 + % Get default optimization parameters. + optimParams = getReweightedOptimizationParams; + end + + if nargin > 3 + % Use ground truth image provided by the user for synthetic data + % experiments. + groundTruth = varargin{1}; + else + % No ground truth available. + groundTruth = []; + end + + if nargout > 2 + % Setup the report structure. + report.SR = {[]}; + report.sigmaNoise = []; + report.sigmaPrior = []; + report.regularizationWeight = []; + report.valError = {[]}; + report.trainError = {[]}; + report.numFunEvals = 0; + report.observationWeights = {}; + report.priorWeights = {}; + end + + %********** Initialization + % Setup the coarse-to-fine optimization parameters. + if ~isempty(optimParams.numCoarseToFineLevels) + % Get user-defined number of coarse-to-fine levels to build the + % image pyramid. + numCoarseToFineLevels = optimParams.numCoarseToFineLevels; + coarseToFineScaleFactors = max([1 (model.magFactor - numCoarseToFineLevels + 1)]) : model.magFactor; + else + % Use default settings for coarse-to-fine optimization. Build image + % pyramid with integer scales between 1 and the desired + % magnification factor. + coarseToFineScaleFactors = min([1 model.magFactor]) : model.magFactor; + end + if isempty(model.SR) + % Initialize super-resolved image by the temporal median of the + % motion-compensated low-resolution frames. + SR = imageToVector( imresize(medfilttemp(LRImages, model.motionParams), coarseToFineScaleFactors(1)) ); + else + % Use the user-defined initial guess. This image needs to be + % resized to the coarsest level of the image pyramid. + SR = imageToVector( imresize(model.SR, coarseToFineScaleFactors(1) * size(LRImages(:,:,1))) ); + end + % Initialize the confidence weights of the observation model. + for frameIdx = 1:size(LRImages, 3) + % Use uniform weights as initial guess. + observationWeights{frameIdx} = ones(numel(LRImages(:,:,frameIdx)), 1); + end + if isempty(model.confidence) + model.confidence = observationWeights; + end + observationWeightsStatic = model.confidence; + + % Iterations for cross validation based hyperparameter selection. + maxCVIter = optimParams.maxCVIter; + + % Decide if the regularization weight should be automatically adjusted + % per iteration. + if isempty(model.imagePrior.weight) + % No user-defined parameter available. Adjust the weight at each + % iteration. + useFixedRegularizationWeight = false; + else + % Use the user-defined regularization weight. + useFixedRegularizationWeight = true; + end + + % Main optimization loop of iteratively re-weighted minimization. + for iter = 1 : optimParams.maxMMIter + + %********** Extract current level from image pyramid. + if iter <= length(coarseToFineScaleFactors) + % Assemble the system matrices for this level of the image + % pyramid. + model.magFactor = coarseToFineScaleFactors(iter); + for frameIdx = 1:size(LRImages, 3) + W{frameIdx} = composeSystemMatrix(size(LRImages(:,:,frameIdx)), model.magFactor, model.psfWidth, model.motionParams{frameIdx}); + Wt{frameIdx} = W{frameIdx}'; + end + + % Propagate the estimate to this level of the image pyramid. + if iter > 1 + SR = imresize(vectorToImage(SR, coarseToFineScaleFactors(iter-1) * size(LRImages(:,:,1))), coarseToFineScaleFactors(iter)*size(LRImages(:,:,1))); + SR = imageToVector(SR); + end + end + % else: we have already reached the desired magnification level. + + %********** Update the observation confidence weights. + sigmaNoise = 0; + if ~isempty(optimParams.observationWeightingFunction) + for frameIdx = 1:size(LRImages, 3) + % Compute the residual error. + y{frameIdx} = imageToVector(LRImages(:,:,frameIdx)); + if ~isempty(model.photometricParams) + if ~isvector(model.photometricParams.mult) + residualError{frameIdx} = getResidualForSingleFrame(SR, y{frameIdx}, W{frameIdx}, model.photometricParams.mult(:,:,frameIdx), model.photometricParams.add(:,:,frameIdx)); + else + residualError{frameIdx} = getResidualForSingleFrame(SR, y{frameIdx}, W{frameIdx}, model.photometricParams.mult(frameIdx), model.photometricParams.add(frameIdx)); + end + else + residualError{frameIdx} = getResidualForSingleFrame(SR, y{frameIdx}, W{frameIdx}); + end + end + [observationWeights, sigmaNoise] = optimParams.observationWeightingFunction.function(residualError, model.confidence, optimParams.observationWeightingFunction.parameters{1:end}); + + % Combine the dynamic observation weights with static, + % user-defined confidence weights. + for frameIdx = 1:size(LRImages, 3) + model.confidence{frameIdx} = observationWeightsStatic{frameIdx} .* observationWeights{frameIdx}; + end + + % else: use uniform weights if no weighting function provided + end + + %********** Update the image prior confidence weights. + sigmaPrior = 0; + if ~isempty(optimParams.priorWeightingFunction) + % Apply sparsity transform to the current estimate of the + % high-resolution image according to the image prior. + model.imagePrior.parameters{1} = model.magFactor * size(LRImages(:,:,1)); + [~, transformedImage] = model.imagePrior.function(SR, model.imagePrior.parameters{1:end-1}); + % Filtering of the sparsity transformed image. + transformedImageFiltered = []; + if iscell(transformedImage) + for l = 1:size(transformedImage,1) + for m = 1:size(transformedImage,2) + if ~isempty(transformedImage{l,m}) + z = vectorToImage( transformedImage{l,m}, model.magFactor * size(LRImages(:,:,1))); + transformedImageFiltered = [transformedImageFiltered; imageToVector( medfilt2(z, [3 3]) )]; + end + end + end + else + transformedImageFiltered = imageToVector( medfilt2(transformedImage, [3 3]) ); + end + if ~exist('priorWeights', 'var') + % Initialize weights at the first iteration. + priorWeights = ones(size(transformedImageFiltered)); + end + if numel(priorWeights) ~= numel(transformedImageFiltered) + % Propagate the weights from the previous level in + % coarse-to-fine optimization to the current level. + priorWeights = imresize(priorWeights, size(transformedImageFiltered)); + end + [priorWeights, sigmaPrior] = optimParams.priorWeightingFunction.function(transformedImageFiltered, priorWeights, optimParams.priorWeightingFunction.parameters{1:end}); + else + % Use uniform weights if no weighting function provided. + priorWeights = 1; + end + model.imagePrior.parameters{end} = priorWeights; + + %********** Hyperparameter selection + if maxCVIter > 1 && ~useFixedRegularizationWeight + % Automatic hyperparameter selection using cross validation. + [model.imagePrior.weight, SR_best] = selectRegularizationWeight; + % Update number of cross validation iterators. + maxCVIter = max([round( 0.5 * maxCVIter ) 1]); + else + % No automatic hyperparameter selection required. + SR_best = SR; + end + + %********** Update estimate for the high-resolution image. + SR_old = SR; + [SR, numFunEvals] = updateHighResolutionImage(SR_best, model, y, W, Wt); + + %********** Check for convergence. + if isConverged(SR, SR_old) && (iter > length(coarseToFineScaleFactors)) + % Convergence tolerance reached. + SR = vectorToImage(SR, model.magFactor * size(LRImages(:,:,1))); + return; + end + + if nargout > 2 + % Log results for current iteration. + report.SR{iter} = vectorToImage(SR, model.magFactor * size(LRImages(:,:,1))); + report.numFunEvals = report.numFunEvals + numFunEvals; + report.sigmaNoise(iter) = sigmaNoise; + report.sigmaPrior(iter) = sigmaPrior; + report.regularizationWeight(iter) = model.imagePrior.weight; + report.observationWeights = cat(1, report.observationWeights, model.confidence); + report.priorWeights = cat(1, report.priorWeights, priorWeights); + if ~isempty(groundTruth) + % Measure PSNR and SSIM for the given ground truth image. + report.psnr(iter) = psnr(imresize(vectorToImage(SR, model.magFactor * size(LRImages(:,:,1))), size(groundTruth)), groundTruth); + report.ssim(iter) = ssim(imresize(vectorToImage(SR, model.magFactor * size(LRImages(:,:,1))), size(groundTruth)), groundTruth); + end + end + + end + + SR = vectorToImage(SR, model.magFactor * size(LRImages(:,:,1))); + + function [SR, numIters] = updateHighResolutionImage(SR, model, y, W, Wt) + + % Setup parameters for SCG optimization. + scgOptions = zeros(1,18); + scgOptions(2) = optimParams.terminationTol; + scgOptions(3) = optimParams.terminationTol; + scgOptions(10) = optimParams.maxSCGIter; + scgOptions(14) = optimParams.maxSCGIter; + + % Perform SCG iterations to update the current estimate of the + % high-resolution image. + if iscolumn(SR) + SR = SR'; + end + [SR, ~, flog] = scg(@imageObjectiveFunc, SR, scgOptions, @imageObjectiveFunc_grad, model, y, W, Wt); + numIters = length(flog); + SR = SR'; + + end + + function [bestLambda, SR_best] = selectRegularizationWeight + + % Split the set of given observations into training and validation + % subset. + for k = 1:size(LRImages, 3) + fractionCvTrainingObservations = optimParams.fractionCVTrainingObservations; + trainObservations{k} = 1 - (randn(size(y{k})) > fractionCvTrainingObservations); %#ok<*AGROW> + y_train{k} = y{k}(trainObservations{k} == 1); + y_val{k} = y{k}(trainObservations{k} == 0); + W_train{k} = W{k}(trainObservations{k} == 1,:); + Wt_train{k} = W_train{k}'; + W_val{k} = W{k}(trainObservations{k} == 0,:); + if ~isempty(model.photometricParams) + if isvector(model.photometricParams.mult) + gamma_m_train{k} = model.photometricParams.mult(k); + gamma_m_val{k} = model.photometricParams.mult(k); + gamma_a_train{k} = model.photometricParams.add(k); + gamma_a_val{k} = model.photometricParams.add(k); + else + gamma_m = model.photometricParams.mult(:,:,k); + gamma_m_train{k} = gamma_m(trainObservations{k} == 1); + gamma_m_val{k} = gamma_m(trainObservations{k} == 0); + gamma_a = model.photometricParams.add(:,:,k); + gamma_a_train{k} = gamma_a(trainObservations{k} == 1); + gamma_a_val{k} = gamma_a(trainObservations{k} == 0); + end + else + gamma_m_train{k} = []; + gamma_m_val{k} = []; + gamma_a_train{k} = []; + gamma_a_val{k} = []; + end + end + + % Setup the model structure for the training subset. + parameterTrainingModel = model; + for k = 1:size(LRImages, 3) + observationConfidenceWeights = model.confidence{k}; + parameterTrainingModel.confidence{k} = observationConfidenceWeights(trainObservations{k} == 1); + end + + % Define search range for adaptive grid search. + if ~isempty(model.imagePrior.weight) + % Refine the search range from the previous iteration. + lambdaSearchRange = logspace(log10(model.imagePrior.weight) - 1/iter, log10(model.imagePrior.weight) + 1/iter, maxCVIter); + else + % Set search range used for initialization. + lambdaSearchRange = logspace(optimParams.hyperparameterCVSearchRange(1), optimParams.hyperparameterCVSearchRange(2), maxCVIter); + bestLambda = median(lambdaSearchRange); + end + + % Perform adaptive grid search over the selected search range. + SR_best = SR; + minValError = Inf; + if exist('report', 'var') + report.valError{iter} = []; + report.trainError{iter} = []; + end + for lambda = lambdaSearchRange + + % Estimate super-resolved image from the training set. + parameterTrainingModel.imagePrior.weight = lambda; + [SR_train, numFunEvals] = updateHighResolutionImage(SR, parameterTrainingModel, y_train, W_train, Wt_train); + + % Determine errors on the training and the validation subset. + valError = 0; + trainError = 0; + for k = 1:size(LRImages, 3) + observationConfidenceWeights = model.confidence{k}; + % Error on the validation subset. + rk_val = getResidualForSingleFrame(SR_train, y_val{k}, W_val{k}, gamma_m_val{k}, gamma_a_val{k}); + valError = valError + sum( observationConfidenceWeights(trainObservations{k} == 0) .* (rk_val.^2) ); + % Error on the training subset. + rk_train = getResidualForSingleFrame(SR_train, y_train{k}, W_train{k}, gamma_m_train{k}, gamma_a_train{k}); + trainError = trainError + sum( observationConfidenceWeights(trainObservations{k} == 1) .* (rk_train.^2) ); + end + if valError < minValError + % Found optimal regularization weight. + bestLambda = lambda; + minValError = valError; + SR_best = SR_train; + end + + if exist('report', 'var') + report.numFunEvals = report.numFunEvals + length(numFunEvals); + % Save errors on training and validation sets. + report.valError{iter} = cat(1, report.valError{iter}, valError); + report.trainError{iter} = cat(1, report.trainError{iter}, trainError); + end + end + + end + + function converged = isConverged(SR, SR_old) + converged = (max(abs(SR_old - SR)) < optimParams.terminationTol); + end + +end + +function f = imageObjectiveFunc(SR, model, y, W, ~) + + if ~iscolumn(SR) + % Reshape to column vector. + SR = SR'; + end + + % Evaluate the data fidelity term. + dataTerm = 0; + for k = 1:length(y) + if ~isempty(model.photometricParams) + if isvector(model.photometricParams.mult) + rk = getResidualForSingleFrame(SR, y{k}, W{k}, model.photometricParams.mult(k), model.photometricParams.add(k)); + else + rk = getResidualForSingleFrame(SR, y{k}, W{k}, model.photometricParams.mult(:,:,k), model.photometricParams.add(:,:,k)); + end + else + rk = getResidualForSingleFrame(SR, y{k}, W{k}); + end + dataTerm = dataTerm + sum( model.confidence{k} .* (rk.^2) ); + end + + % Evaluate image prior for regularization the super-resolved estimate. + priorTerm = model.imagePrior.function(SR, model.imagePrior.parameters{1:end}); + + % Calculate objective function. + f = dataTerm + model.imagePrior.weight * priorTerm; + +end + +function grad = imageObjectiveFunc_grad(SR, model, y, W, Wt) + + if ~iscolumn(SR) + % Reshape to column vector. + SR = SR'; + end + + % Calculate gradient of the data fidelity term w.r.t. the + % super-resolved image. + dataTerm_grad = 0; + for k = 1:length(y) + if ~isempty(model.photometricParams) + if isvector(model.photometricParams.mult) + rk = getResidualForSingleFrame(SR, y{k}, W{k}, model.photometricParams.mult(k), model.photometricParams.add(k)); + else + rk = getResidualForSingleFrame(SR, y{k}, W{k}, model.photometricParams.mult(:,:,k), model.photometricParams.add(:,:,k)); + end + else + rk = getResidualForSingleFrame(SR, y{k}, W{k}); + end + dataTerm_grad = dataTerm_grad - 2*Wt{k} * ( model.confidence{k} .* rk ); + end + + % Calculate gradient of the regularization term w.r.t. the + % super-resolved image. + priorTerm_grad = model.imagePrior.gradient(SR, model.imagePrior.parameters{1:end}); + + % Sum up to total gradient + grad = dataTerm_grad + model.imagePrior.weight * priorTerm_grad; + grad = grad'; + +end + +function r = getResidualForSingleFrame(x, y, W, gamma_m, gamma_a) + + if nargin < 4 || isempty(gamma_m) + gamma_m = 1; + end + if nargin < 5 || isempty(gamma_a) + gamma_a = 0; + end + r = (y - (gamma_m .* (W*x) + gamma_a)); + +end + + \ No newline at end of file diff --git a/teste tensor.py b/teste tensor.py new file mode 100644 index 0000000..44b7991 --- /dev/null +++ b/teste tensor.py @@ -0,0 +1,457 @@ +import SimpleITK as sitk +import os +import numpy as np +from skimage.transform import resize +import json +import numpy as np +from skimage.transform import resize +import subprocess +import scipy.io as sio +from scipy.sparse import csc_matrix, coo_array +# from scipy import sparse +import json +import imresize +from scipy.signal import medfilt2d +from scipy.optimize import fmin_cg + +def smooth_and_resample(image, shrink_factor, smoothing_sigma): + """ + Args: + image: The image we want to resample. + shrink_factor: A number greater than one, such that the new image's size is original_size/shrink_factor. + smoothing_sigma: Sigma for Gaussian smoothing, this is in physical (image spacing) units, not pixels. + Return: + Image which is a result of smoothing the input and then resampling it using the given sigma and shrink factor. + """ + smoothed_image = sitk.SmoothingRecursiveGaussian(image, smoothing_sigma) + + original_spacing = image.GetSpacing() + original_size = image.GetSize() + new_size = [int(sz / float(shrink_factor) + 0.5) for sz in original_size] + new_spacing = [ + ((original_sz - 1) * original_spc) / (new_sz - 1) + for original_sz, original_spc, new_sz in zip( + original_size, original_spacing, new_size + ) + ] + return sitk.Resample( + smoothed_image, + new_size, + sitk.Transform(), + sitk.sitkLinear, + image.GetOrigin(), + new_spacing, + image.GetDirection(), + 0.0, + image.GetPixelID(), + ) + + +def multiscale_demons( + registration_algorithm, + fixed_image, + moving_image, + initial_transform=None, + shrink_factors=None, + smoothing_sigmas=None, +): + """ + Run the given registration algorithm in a multiscale fashion. The original scale should not be given as input as the + original images are implicitly incorporated as the base of the pyramid. + Args: + registration_algorithm: Any registration algorithm that has an Execute(fixed_image, moving_image, displacement_field_image) + method. + fixed_image: Resulting transformation maps points from this image's spatial domain to the moving image spatial domain. + moving_image: Resulting transformation maps points from the fixed_image's spatial domain to this image's spatial domain. + initial_transform: Any SimpleITK transform, used to initialize the displacement field. + shrink_factors: Shrink factors relative to the original image's size. + smoothing_sigmas: Amount of smoothing which is done prior to resmapling the image using the given shrink factor. These + are in physical (image spacing) units. + Returns: + SimpleITK.DisplacementFieldTransform + """ + # Create image pyramid. + fixed_images = [fixed_image] + moving_images = [moving_image] + if shrink_factors: + for shrink_factor, smoothing_sigma in reversed( + list(zip(shrink_factors, smoothing_sigmas)) + ): + fixed_images.append( + smooth_and_resample(fixed_images[0], shrink_factor, smoothing_sigma) + ) + moving_images.append( + smooth_and_resample(moving_images[0], shrink_factor, smoothing_sigma) + ) + + # Create initial displacement field at lowest resolution. + # Currently, the pixel type is required to be sitkVectorFloat64 because of a constraint imposed by the Demons filters. + if initial_transform: + initial_displacement_field = sitk.TransformToDisplacementField( + initial_transform, + sitk.sitkVectorFloat64, + fixed_images[-1].GetSize(), + fixed_images[-1].GetOrigin(), + fixed_images[-1].GetSpacing(), + fixed_images[-1].GetDirection(), + ) + else: + initial_displacement_field = sitk.Image( + fixed_images[-1].GetWidth(), + fixed_images[-1].GetHeight(), + #fixed_images[-1].GetDepth(), + sitk.sitkVectorFloat64, + ) + initial_displacement_field.CopyInformation(fixed_images[-1]) + + # Run the registration. + initial_displacement_field = registration_algorithm.Execute( + fixed_images[-1], moving_images[-1], initial_displacement_field + ) + # Start at the top of the pyramid and work our way down. + for f_image, m_image in reversed( + list(zip(fixed_images[0:-1], moving_images[0:-1])) + ): + initial_displacement_field = sitk.Resample(initial_displacement_field, f_image) + initial_displacement_field = registration_algorithm.Execute( + f_image, m_image, initial_displacement_field + ) + return sitk.DisplacementFieldTransform(initial_displacement_field) + +# Define a simple callback which allows us to monitor registration progress. +def iteration_callback(filter): + print( + "\r{0}: {1:.2f}".format(filter.GetElapsedIterations(), filter.GetMetric()), + end="", + ) + +def imageToVector(img): + return img.flatten('F')[:,np.newaxis] + +def get_residual(SR, LR, W, photometric_params=None): + # GETRESIDUAL Get residual error for low-resolution and super-resolved + # data. + # GETRESIDUAL computes the residual error caused by a super-resolved + # image with the associated low-resolution frames and the model + # parameters (system matrix and photometric parameters). + + if photometric_params is None: + sSR = csc_matrix(SR) + print(LR.shape) + print(W.shape) + print(sSR.shape) + print(type(W)) + print(type(sSR)) + print('não fez') + r = LR - np.multiply(W,sSR) + print('fez') + else: + num_frames = photometric_params.mult.shape[2] + num_lr_pixel = len(LR) // num_frames + + if np.ndim(photometric_params.mult) == 2: # Check if 'mult' is a vector + bm = np.zeros_like(LR) + ba = np.zeros_like(LR) + for k in range(num_frames): + start_idx = (k * num_lr_pixel) + end_idx = ((k + 1) * num_lr_pixel) + bm[start_idx:end_idx] = np.tile(photometric_params.mult[:,:,k], (num_lr_pixel, 1)).flatten() + ba[start_idx:end_idx] = np.tile(photometric_params.mult[:,:,k], (num_lr_pixel, 1)).flatten() + else: + bm = np.zeros_like(LR) + ba = np.zeros_like(LR) + for k in range(num_frames): + start_idx = (k * num_lr_pixel) + end_idx = ((k + 1) * num_lr_pixel) + bm[start_idx:end_idx] = photometric_params.mult[:,:,k].flatten() + ba[start_idx:end_idx] = photometric_params.add[:,:,k].flatten() + + r = LR - bm * np.dot(W, SR) - ba + + return r + +def estimate_observation_confidence_weights(r, weights=None, scale_parameter=None): + # Estimate observation confidence weights + weights_vec = [] + r_vec = [] + r_max = 0.02 # Threshold for median residual to detect outliers + # r_max = 0.1 # Alternative threshold (commented out in MATLAB) + print(len(r)) + print(len(weights)) + + for k in range(len(r)): + r_vec.extend(r[k]) + if np.abs(np.median(r[k])).any() < r_max: + weights_bias_k = 1 + else: + weights_bias_k = 0 + weights_vec.extend([w * weights_bias_k for w in weights[k]]) + + if scale_parameter is None: + if any(w > 0 for w in weights_vec): + # Use adaptive scale parameter estimation + scale_parameter = get_adaptive_scale_parameter(np.array(r_vec)[np.array(weights_vec) > 0], + np.array(weights_vec)[np.array(weights_vec) > 0]) + else: + scale_parameter = 1 + + for k in range(len(r)): + # Estimate local confidence weights for current frame + c = 2 + weights_local = 1 / np.abs(r[k]) + weights_local[np.abs(r[k]) < c * scale_parameter] = 1 / (c * scale_parameter) + weights_local = c * scale_parameter * weights_local + + # Assemble confidence weights from bias weights and local weights + if any(weights_bias_k > 0 for weights_bias_k in weights_vec): + weights[k] = [weights_bias_k * weights_local_i for weights_bias_k, weights_local_i in zip(weights_vec, weights_local)] + else: + weights[k] = [1 * weights_local_i for weights_local_i in weights_local] + + return weights, scale_parameter + +def get_adaptive_scale_parameter(r, weights): + # Nested function to compute adaptive scale parameter + weighted_median_diff = np.abs(r - weighted_median(r, weights)) + scale_parameter = 1.4826 * weighted_median(weighted_median_diff[weights > 0], weights[weights > 0]) + return scale_parameter + +def weighted_median(values, weights): + # Function to compute weighted median + sorted_indices = np.argsort(values) + sorted_values = values[sorted_indices] + sorted_weights = weights[sorted_indices] + cumulative_weights = np.cumsum(sorted_weights) + total_weight = np.sum(sorted_weights) + + median_index = np.searchsorted(cumulative_weights, total_weight / 2.0) + return sorted_values[median_index] + +def get_btv_transformed_image(X, P, alpha0): + # Pad image at the border to perform shift operations + Xpad = np.pad(X, ((P, P), (P, P)), mode='symmetric') + + # Consider shifts in the interval [-P, +P] + btv_transformed_image = [] + rows, cols = X.shape + for l in range(-P, P + 1): + for m in range(-P, P + 1): + if l != 0 or m != 0: + # Shift by l and m pixels + Xshift = Xpad[(P + l):(P + l + rows), (P + m):(P + m + cols)] + + # Apply BTV transformation + transformed_patch = alpha0**(np.abs(l) + np.abs(m)) * (Xshift - X) + + # Apply median filtering + filtered_patch = medfilt2d(transformed_patch, kernel_size=3) + + # Convert the filtered patch to a vector and append to btv_transformed_image + btv_transformed_image.append(filtered_patch.flatten()) + + btv_transformed_image = np.array(btv_transformed_image) + + return btv_transformed_image + +class Parametro: + def __init__(self): + self.P = None + self.alpha = None + self.photometricParams = None + self.imagePrior = None + self.confidence = None + self.SR = None + self.magFactor = None + self.psfWidth = None + self.motionParams = None + self.priori = None + + + def setMotionParams(self,numFrames,registration=False): + if registration == False: + motionParams = [] + for i in range(numFrames): + motionParams.append(np.eye(3)) + self.motionParams = np.array(motionParams) + else: + pass # Registro simultâneo à reconstrução, não implementaremos isso + +model = Parametro() +model.P = 2 +model.alpha = 0.5 +model.magFactor = 2 +model.psfWidth = 0.4 +model.magFactor = 8 +model.motionParams = [[1, 0, 0],[0, 1, 0],[0, 0, 1]] +model.confidence = [] + +class imagePrior(): + def __init__(self): + self.weight = None + +model.imagePrior = imagePrior() +model.imagePrior.weight = 1 + +class optimParams(): + def __init__(self): + self.maxCVIter = None + self.maxMMiter = None + + +SR = [] +confidenceweight = {} +coarseToFineScaleFactors = [*range(min(1,model.magFactor), model.magFactor)] +Y= [] +residualError = [] +W_mat = [] + +folder = 'imagens' + +high_res_images= [] + +for file in os.listdir(folder): + if os.path.isfile(os.path.join(folder,file)): + high_res_images.append(file) + +for img in high_res_images: + + low_res_images= [] + folder_LR = f'LR - {img}' + folder_LR_path = os.path.join(folder,folder_LR) + for lr_file in os.listdir(folder_LR_path): + if os.path.isfile(os.path.join(folder_LR_path,lr_file)): + low_res_images.append(lr_file) + + path_fix = os.path.join(folder_LR_path, low_res_images[0]) + images= [] + images.append(sitk.ReadImage(path_fix, sitk.sitkFloat32)) + + matrices_list = [] + for file_mov in range(1,len(low_res_images)): + path_mov= os.path.join(folder_LR_path, low_res_images[file_mov]) + images.append(sitk.ReadImage(path_mov, sitk.sitkFloat32)) + + fixed_index = 0 + moving_index = file_mov + + fixed_image = images[fixed_index] + moving_image = images[moving_index] + + + # Select a Demons filter and configure it. + demons_filter = sitk.FastSymmetricForcesDemonsRegistrationFilter() + demons_filter.SetNumberOfIterations(20) + # Regularization (update field - viscous, total field - elastic). + demons_filter.SetSmoothDisplacementField(True) + demons_filter.SetStandardDeviations(2.0) + + # Add our simple callback to the registration filter. + demons_filter.AddCommand( + sitk.sitkIterationEvent, lambda: iteration_callback(demons_filter) + ) + + # # Run the registration. + tx = multiscale_demons( + registration_algorithm=demons_filter, + fixed_image=fixed_image, + moving_image=moving_image, + shrink_factors=[1, 1], + smoothing_sigmas=[8, 4], + ) + + moving_resampled = sitk.Resample( + moving_image, + fixed_image, + tx, + sitk.sitkLinear, + 0.0, + moving_image.GetPixelID(), + ) + + matx = sitk.GetArrayFromImage(moving_resampled) + #print(matx) + matrices_list.append(matx) + #############################################################3 + SR = [] + + if len(SR) == 0: + tensor = np.stack(matrices_list,axis=-1) + print(tensor.shape) + sr_med = np.zeros((tensor.shape[0],tensor.shape[1])) + for y in range(tensor.shape[0]): + for x in range(tensor.shape[1]): + a = np.array(tensor[y,x,:]) + for value in a: + if value>0: + np.delete(a,np.where(a==0)) + sr_med[y,x] = np.median(a) + #print(sr_med) + + #Setting Scale + np.array(sr_med) + scale = 1 + img_rsz = resize(sr_med,(sr_med.shape[0],sr_med.shape[1]),anti_aliasing= True ) + #print(img_rsz) + SR = sr_med + + SR = imageToVector(SR) + + for frameIdx in range(0,tensor.shape[2]): + model.confidence.append(np.ones((tensor.shape[0]*tensor.shape[1],1))) + + maxCVIter = optimParams.maxCVIter = 2 + maxMMiter = optimParams.maxMMiter = 2 + + if model.imagePrior.weight is not None: + useFixedRegularizationWeight = True + else: + useFixedRegularizationWeight = False + + for iter in range(0,optimParams.maxMMiter): + W=[] + residualError=[] + if iter <= len(coarseToFineScaleFactors): + model.magFactor = coarseToFineScaleFactors[iter] + + for frameIdx in range(0,tensor.shape[2]): + params = {'imsize': (tensor.shape[0],tensor.shape[1]), 'magFactor': model.magFactor, 'psfWidth' : model.psfWidth, 'motionParams': model.motionParams } + json_object = json.dumps(params, indent=8) + with open("params.json", "w") as outfile: + outfile.write(json_object) + matBin = "/home/labcisne/R2022a/bin/matlab" + # Necessario definir variavel folder se o .py nao estiver na mesma pasta do composeSM.m, CC folder = "" + folder = "/home/labcisne/DocThais/DocThais/Projeto/Testes/SRPython" + # folder = "" + script = "testeComposeSystemMarix.m" + + cmd = [] + cmd.append(matBin) + cmd.append("-sd") + cmd.append(folder) + cmd.append("-batch") + cmd.append(f"run('{script}');") + + subprocess.run(cmd) + W_mat = sio.loadmat('Wmx.mat') + W.append(W_mat['W2'] ) + print('fez 1') + + if iter>0: + print('CHHEEEGGOOOUUUUUUUUUUUUUUUU"') + SR= imresize.imresize(np.reshape(SR,(tensor.shape[0]*coarseToFineScaleFactors[iter-1],tensor.shape[1]*coarseToFineScaleFactors[iter-1]),'C'), + (coarseToFineScaleFactors[iter]/coarseToFineScaleFactors[iter-1])) + SR = imageToVector(SR) + print(f"Elementos lista W: {len(W)}") + for frameIdx in range(0,tensor.shape[2]): + print(f'Imagem número:{frameIdx}\n') + Y.append(imageToVector(tensor[:,:,frameIdx])) + #residualError[frameIdx] = get_residual(SR, Y[frameIdx],W[frameIdx], model.photometricParams) + residualError.append(get_residual(SR, Y[frameIdx],W[frameIdx], model.photometricParams)) + print('fez 2') + + model.confidence, sigmaNoise = estimate_observation_confidence_weights(residualError, model.confidence) + + + + + \ No newline at end of file diff --git "a/teste tradu\303\247\303\243o reweighted.py" "b/teste tradu\303\247\303\243o reweighted.py" index c2a82b0..c643453 100644 --- "a/teste tradu\303\247\303\243o reweighted.py" +++ "b/teste tradu\303\247\303\243o reweighted.py" @@ -6,6 +6,7 @@ import json import imresize from scipy.signal import medfilt2d +from scipy.optimize import fmin_cg # if isempty(model.SR) # % Initialize super-resolved image by the temporal median of the # % motion-compensated low-resolution frames. @@ -35,14 +36,29 @@ # % Use the user-defined regularization weight. # useFixedRegularizationWeight = true; # end -SR = [] -confidenceweight = {} -optimParams = {} -model = {} -coarseToFineScaleFactors = range(min(1,model['magFactor']), model.magFactor) -Y= [] -residualError = [] -W = [] +class Parametro: + def __init__(self): + self.P = None + self.alpha = None + self.photometricParams = None + self.imagePrior = None + self.confidence = None + self.SR = None + self.magFactor = None + self.psfWidth = None + self.motionParams = None + self.priori = None + + + def setMotionParams(self,numFrames,registration=False): + if registration == False: + motionParams = [] + for i in range(numFrames): + motionParams.append(np.eye(3)) + self.motionParams = np.array(motionParams) + else: + pass # Registro simultâneo à reconstrução, não implementaremos isso + def imageToVector(img): return img.flatten('F')[:,np.newaxis] @@ -57,25 +73,25 @@ def get_residual(SR, LR, W, photometric_params=None): if photometric_params is None: r = LR - np.dot(W, SR) else: - num_frames = photometric_params['mult'].shape[2] + num_frames = photometric_params.mult.shape[2] num_lr_pixel = len(LR) // num_frames - if np.ndim(photometric_params['mult']) == 2: # Check if 'mult' is a vector + if np.ndim(photometric_params.mult) == 2: # Check if 'mult' is a vector bm = np.zeros_like(LR) ba = np.zeros_like(LR) for k in range(num_frames): start_idx = (k * num_lr_pixel) end_idx = ((k + 1) * num_lr_pixel) - bm[start_idx:end_idx] = np.tile(photometric_params['mult'][:,:,k], (num_lr_pixel, 1)).flatten() - ba[start_idx:end_idx] = np.tile(photometric_params['add'][:,:,k], (num_lr_pixel, 1)).flatten() + bm[start_idx:end_idx] = np.tile(photometric_params.mult[:,:,k], (num_lr_pixel, 1)).flatten() + ba[start_idx:end_idx] = np.tile(photometric_params.mult[:,:,k], (num_lr_pixel, 1)).flatten() else: bm = np.zeros_like(LR) ba = np.zeros_like(LR) for k in range(num_frames): start_idx = (k * num_lr_pixel) end_idx = ((k + 1) * num_lr_pixel) - bm[start_idx:end_idx] = photometric_params['mult'][:,:,k].flatten() - ba[start_idx:end_idx] = photometric_params['add'][:,:,k].flatten() + bm[start_idx:end_idx] = photometric_params.mult[:,:,k].flatten() + ba[start_idx:end_idx] = photometric_params.add[:,:,k].flatten() r = LR - bm * np.dot(W, SR) - ba @@ -196,6 +212,181 @@ def weighted_median(values, weights): median_index = np.searchsorted(cumulative_weights, total_weight / 2.0) return sorted_values[median_index] +def updateHighResolutionImage(SR, optimParams, model, y, W, Wt): + # Setup parameters for CG optimization. + cgOptions = {'gtol': optimParams['terminationTol'], + 'maxiter': optimParams['maxCGIter'], + 'disp': False} + + # Perform CG iterations to update the current estimate of the + # high-resolution image. + if SR.ndim == 1: + SR = SR[:, np.newaxis] # Convert to column vector if SR is 1D + + def imageObjectiveFunc(x, *args): + # Define your objective function here + # For example: + # return np.sum((y - model @ x)**2) + pass + + def imageObjectiveFunc_grad(x, *args): + # Define your gradient of the objective function here + # For example: + # return -2 * model.T @ (y - model @ x) + pass + + SR, _, _ = fmin_cg(imageObjectiveFunc, SR.flatten(), fprime=imageObjectiveFunc_grad, + args=(model, y, W, Wt), **cgOptions) + + numIters = _['nfev'] # Assuming 'nfev' gives the number of iterations + SR = SR.flatten() + + return SR, numIters + +def isConverged(SR, SR_old, optimParams): + converged = np.max(np.abs(SR_old - SR)) < optimParams['terminationTol'] + return converged + +import numpy as np + +def selectRegularizationWeight(SR, optimParams, model, y, W): + # Initialize variables + maxCVIter = optimParams['maxCVIter'] + fractionCvTrainingObservations = optimParams['fractionCVTrainingObservations'] + hyperparameterCVSearchRange = optimParams['hyperparameterCVSearchRange'] + + bestLambda = None + SR_best = SR + minValError = np.inf + + # Split the set of given observations into training and validation subset. + trainObservations = [] + y_train = [] + y_val = [] + W_train = [] + Wt_train = [] + W_val = [] + + for k in range(len(y)): + trainObservations_k = np.random.rand(len(y[k])) > fractionCvTrainingObservations + trainObservations.append(trainObservations_k) + y_train.append(y[k][trainObservations_k]) + y_val.append(y[k][~trainObservations_k]) + W_train.append(W[k][trainObservations_k, :]) + Wt_train.append(W_train[k].T) + W_val.append(W[k][~trainObservations_k, :]) + + # Setup the model structure for the training subset. + parameterTrainingModel = model.copy() + for k in range(len(y)): + parameterTrainingModel['confidence'][k] = model['confidence'][k][trainObservations[k]] + + # Define search range for adaptive grid search. + if 'imagePrior' in model and 'weight' in model['imagePrior'] and model['imagePrior']['weight'] is not None: + # Refine the search range from the previous iteration. + lambdaSearchRange = np.logspace(np.log10(model['imagePrior']['weight']) - 1/maxCVIter, + np.log10(model['imagePrior']['weight']) + 1/maxCVIter, + maxCVIter) + else: + # Set search range used for initialization. + lambdaSearchRange = np.logspace(hyperparameterCVSearchRange[0], hyperparameterCVSearchRange[1], maxCVIter) + bestLambda = np.median(lambdaSearchRange) + + # Perform adaptive grid search over the selected search range. + for lambda_ in lambdaSearchRange: + # Estimate super-resolved image from the training set. + parameterTrainingModel['imagePrior']['weight'] = lambda_ + SR_train, numFunEvals = updateHighResolutionImage(SR, parameterTrainingModel, y_train, W_train, Wt_train) + + # Determine errors on the training and the validation subset. + valError = 0 + trainError = 0 + for k in range(len(y)): + observationConfidenceWeights = model['confidence'][k] + # Error on the validation subset. + valError += np.sum(observationConfidenceWeights[~trainObservations[k]] * (W_val[k] @ SR_train - y_val[k])**2) + # Error on the training subset. + trainError += np.sum(observationConfidenceWeights[trainObservations[k]] * (W_train[k] @ SR_train - y_train[k])**2) + + if valError < minValError: + # Found optimal regularization weight. + bestLambda = lambda_ + minValError = valError + SR_best = SR_train + + # Save errors on training and validation sets if requested + # (assuming report is a dictionary where these values are stored) + + + return bestLambda, SR_best + +def imageObjectiveFunc(SR, model, y, W, _): + import numpy as np + + if not np.ndim(SR) == 1: + # Reshape to column vector. + SR = SR.reshape(-1, 1) + + # Evaluate the data fidelity term. + dataTerm = 0 + for k in range(len(y)): + dataTerm += np.sum(model['confidence'][k] * (y[k] - W[k] @ SR)**2) + + # Evaluate image prior for regularization the super-resolved estimate. + priorTerm = model['imagePrior']['function'](SR, *model['imagePrior']['parameters']) + + # Calculate objective function. + f = dataTerm + model['imagePrior']['weight'] * priorTerm + + return f + +def imageObjectiveFunc_grad(SR, model, y, W, Wt): + import numpy as np + + if not np.ndim(SR) == 1: + # Reshape to column vector. + SR = SR.reshape(-1, 1) + + # Calculate gradient of the data fidelity term w.r.t. SR. + dataTerm_grad = 0 + for k in range(len(y)): + dataTerm_grad -= 2 * Wt[k] @ (model['confidence'][k] * (y[k] - W[k] @ SR)) + + # Calculate gradient of the regularization term w.r.t. SR. + priorTerm_grad = model['imagePrior']['gradient'](SR, *model['imagePrior']['parameters']) + + # Sum up to total gradient + grad = dataTerm_grad + model['imagePrior']['weight'] * priorTerm_grad + grad = grad.flatten() # Ensure gradient is returned as a 1D array + + return grad + + +model = Parametro() +model.P = 2 +model.alpha = 0.5 +model.magFactor = 2 +model.psfWidth = 0.4 + +class imagePrior(): + def __init__(self): + self.weight = None + +model.imagePrior = imagePrior() +model.imagePrior.weight = 1 + +class optimParams(): + def __init__(self): + self.maxCVIter = None + self.maxMMiter = None + +SR = [] +confidenceweight = {} +coarseToFineScaleFactors = range(min(1,model.magFactor), model.magFactor) +Y= [] +residualError = [] +W = [] + if len(SR) == 0: tensor = np.stack(matrices_list,axis=-1) @@ -222,24 +413,23 @@ def weighted_median(values, weights): for frameIdx in range(0,tensor.shape[2]): confidenceweight[frameIdx] = np.ones(len(tensor.shape[:,:,frameIdx],1)) -maxCVIter = optimParams['maxCVIter'] +maxCVIter = optimParams.maxCVIter + -if model['imagePrior'] and model['imagePrior'].get('weight', None) is not None: +if model.imagePrior.weight is not None: useFixedRegularizationWeight = True else: useFixedRegularizationWeight = False -for iter in range(0,optimParams['maxMMiter']): +for iter in range(0,optimParams.maxMMiter): if iter <= len(coarseToFineScaleFactors): - model['magFactor'] = coarseToFineScaleFactors[iter] + model.magFactor = coarseToFineScaleFactors[iter] for frameIdx in range(0,tensor.shape[2]): - params = {'size': tensor.shape[:,:,frameIdx], 'magFactor': model['magFactor'], 'psfWidth' : model['psfWidth'], 'motionParams': model['motionParams'] } + params = {'size': tensor.shape[:,:,frameIdx], 'magFactor': model.magFactor, 'psfWidth' : model.psfWidth, 'motionParams': model.motionParams } json_object = json.dumps(params, indent=8) with open("params.json", "w") as outfile: outfile.write(json_object) - - matBin = "/home/labcisne/R2022a/bin/matlab" # Necessario definir variavel folder se o .py nao estiver na mesma pasta do composeSM.m, CC folder = "" folder = "/home/labcisne/DocThais/DocThais/Projeto/Testes/SRPython" @@ -291,4 +481,7 @@ def weighted_median(values, weights): SR_best = SR SR_old = SR - SR, numFumEvals = \ No newline at end of file + SR, numFumEvals = updateHighResolutionImage(SR_best, model, y, W, Wt) + + if isConverged(SR, SR_old) and iter>len(coarseToFineScaleFactors): + np.reshape(SR,(tensor.shape[0]*model['magFactor'],tensor.shape[1]*model['magFactor']),'C')