diff --git a/AMaSiNe User Manual 07272020.pdf b/AMaSiNe User Manual 07272020.pdf new file mode 100644 index 0000000..9802da8 Binary files /dev/null and b/AMaSiNe User Manual 07272020.pdf differ diff --git a/STEP_0_Parameters.m b/STEP_0_Parameters.m new file mode 100644 index 0000000..93b0578 --- /dev/null +++ b/STEP_0_Parameters.m @@ -0,0 +1,39 @@ +%% 1. Image Directory +main_folder_dir='C:\Users\junhs\Downloads\AMaSiNe-master\AMaSiNe-with_Elastix_07162020'; +cd(main_folder_dir); +addpath(genpath(main_folder_dir)); + +%% 2. Image names and slice order +img_format ='jpg'; % if tif = 'tif' ; jpg='jpg' +Slice_AP_orPA= 1; % If the brain is sliced from anterior to posterior, set this value = 1 + % posterior to anterior, set this value = -1 +slide_digit=3; +scene_digit=4; +channel_digit= 5; + +%% 3. Anchor Image IDs for Angle Finding (for STEP_2 and 3) +anc_img_IDs= sort([1 4]); +img_IDs_reBoundary=[2]; +threshold_scale = 1.5; %for step 3 only - increase this number if your slice boundary is smaller than what you expect; + %decrease this number if slice + %boundary is not clean enough. + +%% 4. Image Parameters +xy_pix=0.653 * 2; % Pixel size = um/pixel +Name_Channels={'eGFP','DAPI'}; % In the right order +Color_Channel_Structure=2; % ID of Color Channel to be used for angle finding process + % DAPI or NISSL is very strongly recommended for + % the angle finding process (Step_4) +Structure_stain={'DAPI'}; % Choose one of the three : 'DAPI','AutoF','Nissl' +Color_Channel_Interest=[1]; % Color channel in which labelled cells are imaged (e.g. eGFP, tdTomato) + +%% 5. Detection Parameters +soma_radius=[10 16]; % range of "radius" of labelled soma (in um) to be searched for + +cell_det_thresh = 0.25; % Intensity difference between a cell and its background for a cell to be detected as a cell + % Lower this value, you'd get a better chance + % of detecting cells dim, but you also risk detecting noise as a cell + +%% 6. Allen Atlas Info +size_vol = [528 320 456]; % Matrix size of the Allen 3D atlas +ref_atlas_vox_res=25; % Allen 3D Ref atlas : 1 voxel= 25um x 25um x 25um diff --git a/STEP_2_Anchor_Image_Selection.m b/STEP_2_Anchor_Image_Selection.m index 1b2ee73..b4aed84 100644 --- a/STEP_2_Anchor_Image_Selection.m +++ b/STEP_2_Anchor_Image_Selection.m @@ -1,8 +1,8 @@ warning('off') -STEP_0_Parameters; +addpath(genpath(pwd)) img_name=Img_filename_list; - load('Step1_Outline_result.mat') +STEP_0_Parameters; ref_resc=1/0.5; downscaled_xy_pix=ref_atlas_vox_res*ref_resc; diff --git a/STEP_3_Optional_Redraw_Slice_Boundary.m b/STEP_3_Optional_Redraw_Slice_Boundary.m index 962362f..f8ba180 100644 --- a/STEP_3_Optional_Redraw_Slice_Boundary.m +++ b/STEP_3_Optional_Redraw_Slice_Boundary.m @@ -1,4 +1,5 @@ clear all; clc; +addpath(genpath(pwd)) load('Step1_Outline_result.mat') warning('off') STEP_0_Parameters; diff --git a/STEP_4_Angle_Finder.m b/STEP_4_Angle_Finder.m index 9c9d64d..7f2590b 100644 --- a/STEP_4_Angle_Finder.m +++ b/STEP_4_Angle_Finder.m @@ -1,3 +1,4 @@ +clear all; close all; warning('off') STEP_0_Parameters; img_name=Img_filename_list; diff --git a/STEP_4to5_Cell_Detection.m b/STEP_4to5_Cell_Detection.m index 371d155..149ecf7 100644 --- a/STEP_4to5_Cell_Detection.m +++ b/STEP_4to5_Cell_Detection.m @@ -1,14 +1,13 @@ profile on -close all -clear all +close all; clear all +addpath(genpath(pwd)); warning('off') -STEP_0_Parameters; load('ANO_roi_edge'); load('Step_4_Angle_Search_Result'); +STEP_0_Parameters; img_name=Img_filename_list; -mkdir Image_Analysed -mkdir Image_Analysed_ROI_absent +mkdir Image_Analysed_ROI_absent_before_warping %% Prepartion : indices for both the anchor and non-anchor imgs img_idx=anc_img_IDs; @@ -18,7 +17,8 @@ img_AP = ap_found; else for img_ID=1:length(img_idx)-1 - img_AP=[img_AP, linspace(ap_found(img_ID),ap_found(img_ID+1),img_idx(img_ID+1)-img_idx(img_ID)+1)]; + img_AP=[img_AP, linspace(ap_found(img_ID),... + ap_found(img_ID+1),img_idx(img_ID+1)-img_idx(img_ID)+1)]; end end @@ -28,42 +28,46 @@ else img_idx=max(img_idx):-1:min(img_idx); end -xy_pix_resc_factor = 1; + %% MAIN PART -for img_ID=1:length(img_AP) - disp(['Cell Detection: ', num2str(img_ID), '/', num2str(length(img_AP))]); +parfor img_ID=1:length(img_AP) + disp(['Cell Detection: ', num2str(img_ID), '/', ... + num2str(length(img_AP))]); img_ID for color_ch_ID=1:length(Color_Channel_Interest) if ~isempty(img_info(img_idx(img_ID)).slice_window) cell_detection_rs(img_ID).img_AP_pos= img_AP(img_ID); - img_Color=(imread(img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)})); + img_Color=(imread(img_name{img_idx(img_ID),... + Color_Channel_Interest(color_ch_ID)})); try img_Color=rgb2gray(img_Color); end img_Color_pad=zeros(size(img_Color)); - img_Color=img_Color(img_info(img_idx(img_ID)).slice_window(1):img_info(img_idx(img_ID)).slice_window(2),... - img_info(img_idx(img_ID)).slice_window(3):img_info(img_idx(img_ID)).slice_window(4)); - - img_Color_pad=img_Color_pad(img_info(img_idx(img_ID)).slice_window(1):img_info(img_idx(img_ID)).slice_window(2),... - img_info(img_idx(img_ID)).slice_window(3):img_info(img_idx(img_ID)).slice_window(4)); + img_Color=img_Color(img_info(img_idx(img_ID)).slice_window(1):... + img_info(img_idx(img_ID)).slice_window(2),... + img_info(img_idx(img_ID)).slice_window(3):... + img_info(img_idx(img_ID)).slice_window(4)); + img_Color_pad=img_Color_pad(img_info(img_idx(img_ID)).slice_window(1):... + img_info(img_idx(img_ID)).slice_window(2),... + img_info(img_idx(img_ID)).slice_window(3):... + img_info(img_idx(img_ID)).slice_window(4)); img_Color_pad(img_info(img_idx(img_ID)).bnd_pix_ind)=1; img_Color_pad=(imfill(img_Color_pad)); img_Color_pad=uint8(logical(img_Color_pad)); img_Color=img_Color.*img_Color_pad; img_Color=padarray(img_Color,round([3000 3000]/(xy_pix))); - - atlas_exp_scale_fac=(25/xy_pix)*xy_pix_resc_factor; - - img_Color=imresize(img_Color,xy_pix_resc_factor); - + %%% detect cells across the whole slice image - cell_detected_all= SomaDetection0827 (img_Color, xy_pix, soma_radius,cell_det_thresh); + cell_detected_all= SomaDetection0827(... + img_Color, xy_pix, soma_radius,cell_det_thresh); if ~isempty(cell_detected_all) - out_bnd_alpha = ref_boundarypad_0809_step5( img_Color, xy_pix/xy_pix_resc_factor ); - out_bnd=inShape(out_bnd_alpha,cell_detected_all(:,2),cell_detected_all(:,1)); + out_bnd_alpha = ref_boundarypad_0809_step5(... + img_Color, xy_pix); + out_bnd=inShape(out_bnd_alpha,cell_detected_all(:,2),... + cell_detected_all(:,1)); cell_detected_all=cell_detected_all(~out_bnd,:); end @@ -76,7 +80,8 @@ cell_detected_all_pos_ind=[]; end - cell_detection_rs(img_ID).Color_Cells(color_ch_ID).cell_locations=cell_detected_all; + cell_detection_rs(img_ID).Color_Cells(color_ch_ID).cell_locations=... + cell_detected_all; figure; imshow(img_Color); hold on; @@ -84,19 +89,21 @@ image_analyzed=figure; imshow(img_Color,[]); hold on - title(strcat({'Image Name : '}, img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},... - {' - Ch No. '},num2str(color_ch_ID)),'Interpreter', 'none'); + title(strcat({'Image Name : '}, img_name{img_idx(img_ID),... + Color_Channel_Interest(color_ch_ID)},... + {' - Ch No. '},num2str(color_ch_ID)),... + 'Interpreter', 'none'); if ~isempty(cell_detected_all) - scatter(cell_detected_all(:,1),cell_detected_all(:,2),9,'r','filled') + scatter(cell_detected_all(:,1),... + cell_detected_all(:,2),9,'r','filled') end - save_name=strcat('/Image_Analysed_ROI_absent/',img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},'.fig'); + save_name=strcat('/Image_Analysed_ROI_absent_before_warping/',... + img_name{img_idx(img_ID),... + Color_Channel_Interest(color_ch_ID)},'.fig'); saveas( image_analyzed,[pwd save_name]); close(image_analyzed) -% % % [B,RB]=imwarp(A,tform); -% % % [xdataT,ydataT]=transformPointsForward(tform,xdata,ydata); -% % % [xdataI,ydataI]=worldToIntrinsic(RB,xdataT,ydataT); end end end diff --git a/STEP_5_Transform_and_ROI_drawing.m b/STEP_5_Transform_and_ROI_drawing.m index 436edd8..ac2ae87 100644 --- a/STEP_5_Transform_and_ROI_drawing.m +++ b/STEP_5_Transform_and_ROI_drawing.m @@ -1,15 +1,21 @@ profile on -close all -clear all +close all; clear all +addpath(genpath(pwd)); warning('off') -STEP_0_Parameters; load('ANO_roi_edge'); load('Step_4_Angle_Search_Result'); +load('Step_4to5_Cell_Detection_Result.mat') %% Load Cell detection From STEP4to5 +STEP_0_Parameters; + img_name=Img_filename_list; mkdir Image_Analysed mkdir Image_Analysed_ROI_absent +elastixParams = { + strcat(pwd,'\Core_Functions\elastix_functions\warping_parameters_Affine.txt'),... + strcat(pwd,'\Core_Functions\elastix_functions\warping_parameters_BSpline.txt')}; %%% elastix + %% Prepartion : indices for both the anchor and non-anchor imgs img_idx=anc_img_IDs; ap_found=max_APpos_stage_final; @@ -18,7 +24,8 @@ img_AP = ap_found; else for img_ID=1:length(img_idx)-1 - img_AP=[img_AP, linspace(ap_found(img_ID),ap_found(img_ID+1),img_idx(img_ID+1)-img_idx(img_ID)+1)]; + img_AP=[img_AP, linspace(ap_found(img_ID),ap_found(img_ID+1),... + img_idx(img_ID+1)-img_idx(img_ID)+1)]; end end @@ -50,33 +57,37 @@ VOL_rot=imwarp(VOL,tf_atlas,'cubic'); ANO_rot=imwarp(ANO_roi_edge,tf_atlas,'nearest'); -%% Load Cell detection From STEP4to5 -load('Step_4to5_Cell_Detection_Result.mat') + %% MAIN PART parfor_progress(length(img_AP)); errorneous=false(1,length(img_AP)); +img_essence(1:length(img_AP))=struct; -for img_ID=1:length(img_AP) +parfor img_ID=1:length(img_AP) img_ID if ~isempty(img_info(img_idx(img_ID)).slice_window) img_essence(img_ID).img_AP_pos= img_AP(img_ID); current_ap=img_AP(img_ID)-size(VOL,3)/2; - current_ap=round(size(VOL_rot,3)/2+current_ap*cosd(pitch_found)*cosd(yaw_found)); + current_ap=round(size(VOL_rot,3)/2+current_ap*... + cosd(pitch_found)*cosd(yaw_found)); %%%%%%%%%%%%%%%%%%%%%%%% Atlas Slice Prepartion for Transformation %%%%%%%%%%%%%%%%%%%% tform_general_resc_factor=1; img_ref=uint8(squeeze(VOL_rot(:,:,current_ap))); img_ref=imadjust(img_ref,stretchlim(img_ref,0),[0 1]); img_ref=imresize(img_ref,tform_general_resc_factor); - img_ref=padarray(img_ref,round([3000 3000]/(ref_atlas_vox_res/tform_general_resc_factor))); + img_ref=padarray(img_ref,round([3000 3000]/... + (ref_atlas_vox_res/tform_general_resc_factor))); %%%%%%%%%%%%%%%%%%%%%%%%%%% Get Annotated Slice %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% img_ANO=uint8(squeeze(ANO_rot(:,:,current_ap))); img_ANO=padarray(img_ANO,round([3000 3000]/(ref_atlas_vox_res))); - img_ANO=imresize(img_ANO,xy_pix_resc_factor*(ref_atlas_vox_res/xy_pix),'method','nearest','Antialiasing',false); + img_ANO_downscale=img_ANO; + img_ANO=imresize(img_ANO,(ref_atlas_vox_res/xy_pix),... + 'method','nearest','Antialiasing',false); %%%%%%%%%%%%%%%%%%%%%%%% Slice Preparation for Transformation %%%%%%%%%%%%%%%%%%%%%%%%%% @@ -86,333 +97,215 @@ end img_act_pad=zeros(size(img_act)); - img_act=img_act(img_info(img_idx(img_ID)).slice_window(1):img_info(img_idx(img_ID)).slice_window(2),... - img_info(img_idx(img_ID)).slice_window(3):img_info(img_idx(img_ID)).slice_window(4)); - - img_act_pad=img_act_pad(img_info(img_idx(img_ID)).slice_window(1):img_info(img_idx(img_ID)).slice_window(2),... - img_info(img_idx(img_ID)).slice_window(3):img_info(img_idx(img_ID)).slice_window(4)); + img_act=img_act(img_info(img_idx(img_ID)).slice_window(1):... + img_info(img_idx(img_ID)).slice_window(2),... + img_info(img_idx(img_ID)).slice_window(3):... + img_info(img_idx(img_ID)).slice_window(4)); + img_act_pad=img_act_pad(img_info(img_idx(img_ID)).slice_window(1):... + img_info(img_idx(img_ID)).slice_window(2),... + img_info(img_idx(img_ID)).slice_window(3):... + img_info(img_idx(img_ID)).slice_window(4)); img_act_pad(img_info(img_idx(img_ID)).bnd_pix_ind)=1; img_act_pad=(imfill(img_act_pad)); img_act_pad=uint8(logical(img_act_pad)); - img_act=img_act.*(img_act_pad); img_act=padarray(img_act,round([3000 3000]/(xy_pix))); - img_act_reserve=img_act; - img_act=imadjust(img_act,stretchlim(img_act,0.00),[0 1]); + img_act=img_act+( img_act-imgaussfilt(img_act,0.5*201,... + 'FilterSize',[3 3]*603,'FilterDomain','frequency'))*5; + img_act=imresize(img_act,xy_pix/ref_atlas_vox_res*... + tform_general_resc_factor); + - img_act=img_act+( img_act-imgaussfilt( img_act,0.5*201,'FilterSize',[3 3]*603,... - 'FilterDomain','frequency'))*5; + %%%%%%%%%%%%%%%%%%%%%%%%%%% Find Transformation Matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + elastixDir_currImg = fullfile(pwd,... + strcat('\img_warping_log\img_ID_',num2str(img_ID))); %%% elastix + mkdir(elastixDir_currImg); - img_act=imresize(img_act,xy_pix/ref_atlas_vox_res*tform_general_resc_factor); - %%%%%%%%%%%%%%%%%%%%%%%%%%% Find Transformation Matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% try - [matched_img, tform_1, tform_2, tform_3, ... - matchedPoints_act_400, matchedPoints_ref_400, ... - matchedPoints_act_200, matchedPoints_ref_200, matchedPoints_act_100, matchedPoints_ref_100, ... - matchedPoints_act_50, matchedPoints_ref_50] = ... - tform_finder1015( img_act, img_ref, ref_atlas_vox_res/tform_general_resc_factor, [13 13] ); - catch - errorneous(img_ID)=true; - continue; + [matched_img,transform_params]=elastix(img_act, img_ref,elastixDir_currImg,elastixParams); + matched_img=uint8(matched_img); end - if isnan(matchedPoints_act_50) %% if image matching is unsuccessful, skip current iteration - disp(strcat({'Not enough matching points for img #'},num2str(img_ID))) + if isempty(matched_img) || isempty(transform_params) + errorneous(img_ID)=true; continue; end - T1=tform_1.T; T2=tform_2.T; T3=tform_3.T; - tform_combined=mtimes(T1, T2); - tform_combined=mtimes(tform_combined, T3); - tf_general = affine2d(tform_combined); + img_essence(img_ID).transform_params_downscaled=transform_params; - %%%%%%%%%%%%%%%%%%%%%%%% Find translation %%%%%%%%%%%%%%%%%%%%%%%%% - img_act_trans=imwarp(img_act,tf_general); - trans_xcorr = conv2(matched_img,rot90(conj(img_act_trans),2)); - [max_cc, imax] = max(abs(trans_xcorr(:))); - [ypeak, xpeak] = ind2sub(size(trans_xcorr),imax(1)); - translation_vec=([ypeak, xpeak]-size(img_act_trans)); - translation_vec=[translation_vec(2), translation_vec(1)]; - translation_vec=translation_vec*(ref_atlas_vox_res/tform_general_resc_factor)/xy_pix; + %%%%%%%%%%%%%%%%%%%%%%%%% Get DownScaled Img for 3D Rec. %%%%%%%%%%%%%%%%%%%%%%% - img_essence(img_ID).tform.affine=tf_general; - img_essence(img_ID).tform.trans=translation_vec; + img_act_downscale = matched_img; + %%% Match size of images + whos_bigger=size(img_ANO_downscale)-size(img_act_downscale); + if whos_bigger(1)>=0; x_dim=size(img_ANO_downscale,1); + else; x_dim=size(img_act_downscale,1); end + if whos_bigger(2)>=0; y_dim=size(img_ANO_downscale,2); + else; y_dim=size(img_act_downscale,2); end - %%%%%%%%%%%%%%%%%%%%%%%%% Get DownScaled Img for 3D Rec. %%%%%%%%%%%%%%%%%%%%%%% - img_downscaled=(imread(img_name{img_idx(img_ID), Color_Channel_Structure})); - try - img_downscaled=rgb2gray(img_downscaled); - end - - img_downscaled=img_downscaled(img_info(img_idx(img_ID)).slice_window(1):img_info(img_idx(img_ID)).slice_window(2),... - img_info(img_idx(img_ID)).slice_window(3):img_info(img_idx(img_ID)).slice_window(4)); - img_downscaled=img_downscaled.*img_act_pad; - img_downscaled=padarray(img_downscaled,round([3000 3000]/(xy_pix))); + img_act_sizeMatch=[x_dim,y_dim]-size(img_act_downscale); + img_act_sizeMatch(img_act_sizeMatch<0)=0; + img_ANO_sizeMatch=[x_dim,y_dim]-size(img_ANO_downscale); + img_ANO_sizeMatch(img_ANO_sizeMatch<0)=0; + img_act_downscale=padarray(img_act_downscale,... + img_act_sizeMatch,'post'); - img_act_downscale = imwarp(img_downscaled,tf_general); - img_act_downscale = imtranslate(img_act_downscale,translation_vec); + img_essence(img_ID).transformed_img_downscaled=... + img_act_downscale; - %%% Match size of images - whos_bigger=size(img_ANO)-size(img_act_downscale); - if whos_bigger(1)>=0; x_dim=size(img_ANO,1); else; x_dim=size(img_act_downscale,1); end - if whos_bigger(2)>=0; y_dim=size(img_ANO,2); else; y_dim=size(img_act_downscale,2); end + %%%%%%%%%%%% Warp Images in Their Original Pix Resolution %%%%%%%%%%%%%%%%%%%%%%% - img_act_sizeMatch=[x_dim,y_dim]-size(img_act_downscale); img_act_sizeMatch(img_act_sizeMatch<0)=0; - img_ANO_sizeMatch=[x_dim,y_dim]-size(img_ANO); img_ANO_sizeMatch(img_ANO_sizeMatch<0)=0; + %%% modify transform parameters for warping IMAGES %%% + mov_scale_factor=(ref_atlas_vox_res/xy_pix); - img_act_downscale=padarray(img_act_downscale,img_act_sizeMatch,'post'); + transform_params_img_ori_scale=... + img_essence(img_ID).transform_params_downscaled + transform_params_img_ori_scale.TransformParameters{1, 1}.Size=... + fliplr(size(img_ANO)); %% elastix vs matlab... + transform_params_img_ori_scale.TransformParameters{1, 1}.TransformParameters(5:6) = ... + transform_params_img_ori_scale.TransformParameters{1, 1}.TransformParameters(5:6)*mov_scale_factor; + transform_params_img_ori_scale.TransformParameters{1, 1}.CenterOfRotationPoint = ... + round((fliplr(size(img_ref))/2)*mov_scale_factor); + transform_params_img_ori_scale.TransformParameters{1, 2}.Size=... + fliplr(size(img_ANO)); %% elastix vs matlab... + transform_params_img_ori_scale.TransformParameters{1, 2}.GridSpacing = ... + transform_params_img_ori_scale.TransformParameters{1, 2}.GridSpacing*mov_scale_factor; + transform_params_img_ori_scale.TransformParameters{1, 2}.GridOrigin = ... + transform_params_img_ori_scale.TransformParameters{1, 2}.GridOrigin*mov_scale_factor; + transform_params_img_ori_scale.TransformParameters{1, 2}.TransformParameters = ... + transform_params_img_ori_scale.TransformParameters{1, 2}.TransformParameters*mov_scale_factor; - img_essence(img_ID).image_downscaled=imresize(img_act_downscale,xy_pix/ref_atlas_vox_res); - - %%%%%%%%%%%%%%%%%%%%%%%%% Find Cells in ROI %%%%%%%%%%%%%%%%%%%%%%% - if max(max(img_ANO))~=0 - %%% Find Cells in Each ROI %%% - for color_ch_ID=1:length(Color_Channel_Interest) - cells_in_this_image = cell_detection_rs(img_ID).Color_Cells(color_ch_ID).cell_locations; - img_Color=(imread(img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)})); - try - img_Color=rgb2gray(img_Color); - end - - img_Color=img_Color(img_info(img_idx(img_ID)).slice_window(1):img_info(img_idx(img_ID)).slice_window(2),... - img_info(img_idx(img_ID)).slice_window(3):img_info(img_idx(img_ID)).slice_window(4)); - img_Color=img_Color.*img_act_pad; - img_Color=padarray(img_Color,round([3000 3000]/(xy_pix))); - - [img_act_no_scale, tf_general_RB] = imwarp(img_Color,tf_general); - [xdataT_cell, ydataT_cell] = transformPointsForward(tf_general, cells_in_this_image(:,1), cells_in_this_image(:,2)); - [xdataI_cell, ydataI_cell] = worldToIntrinsic(tf_general_RB, xdataT_cell, ydataT_cell); - - [img_act_no_scale, tf_translation_RB] = imtranslate(img_act_no_scale,translation_vec); - cell_warped_x = xdataI_cell+translation_vec(1); - cell_warped_y = ydataI_cell+translation_vec(2); - - resc_warp = size(img_act_no_scale); - - atlas_exp_scale_fac=(25/xy_pix)*xy_pix_resc_factor; - - img_act_no_scale=imresize(img_act_no_scale,xy_pix_resc_factor); - cell_warped_x = cell_warped_x * xy_pix_resc_factor; - cell_warped_y = cell_warped_y * xy_pix_resc_factor; - - img_ref_no_scale=imresize(img_ref,atlas_exp_scale_fac); - - outputView_no_scale = imref2d(size(img_ref_no_scale)); - - matchedPoints_act_400_tuned=matchedPoints_act_400*atlas_exp_scale_fac; - matchedPoints_ref_400_tuned=matchedPoints_ref_400*atlas_exp_scale_fac; - - matchedPoints_act_200_tuned=matchedPoints_act_200*atlas_exp_scale_fac; - matchedPoints_ref_200_tuned=matchedPoints_ref_200*atlas_exp_scale_fac; - - matchedPoints_act_100_tuned=matchedPoints_act_100*atlas_exp_scale_fac; - matchedPoints_ref_100_tuned=matchedPoints_ref_100*atlas_exp_scale_fac; - - matchedPoints_act_50_tuned=matchedPoints_act_50*atlas_exp_scale_fac; - matchedPoints_ref_50_tuned=matchedPoints_ref_50*atlas_exp_scale_fac; - - try - tform_noScale_400 = fitgeotrans(matchedPoints_act_400_tuned,matchedPoints_ref_400_tuned,'pwl'); - [img_warped_no_scale,~] = imwarp_custom((img_act_no_scale), tform_noScale_400,... - 'cubic','OutputView',outputView_no_scale,'WarpRescale',resc_warp); - if ~isempty(cell_warped_x) - distfun = @(p) norm(transformPointsInverse(tform_noScale_400, p) - double([cell_warped_x, cell_warped_y])); - options = optimset('Display','none','MaxFunEvals', numel(cell_warped_x)*500, 'MaxIter', numel(cell_warped_x)*500); - [cell_warpedXY, err] = fminsearch(distfun, double([cell_warped_x, cell_warped_y]), options); - cell_warped_x = cell_warpedXY(:,1); - cell_warped_y = cell_warpedXY(:,2); - end - catch -% disp([img_name{img_idx(img_ID), Color_Channel_Structure}, ': PWL fails in 400um scale - Trying LWM']); -% try -% tform_noScale_400 = fitgeotrans(matchedPoints_act_400_tuned,matchedPoints_ref_400_tuned,'lwm', 12); -% [ img_warped_no_scale,~] = imwarp_custom((img_act_no_scale), tform_noScale_400,... -% 'cubic','OutputView',outputView_no_scale,'WarpRescale',resc_warp); -% if ~isempty(cell_warped_x) -% distfun = @(p) norm(transformPointsInverse(tform_noScale_400, p) - double([cell_warped_x, cell_warped_y])); -% options = optimset('Display','none','MaxFunEvals', numel(cell_warped_x)*500, 'MaxIter', numel(cell_warped_x)*500); -% [cell_warpedXY, err] = fminsearch(distfun, double([cell_warped_x, cell_warped_y]), options); -% cell_warped_x = cell_warpedXY(:,1); -% cell_warped_y = cell_warpedXY(:,2); -% end -% catch - disp(['Error during analyzing : ', img_name{img_idx(img_ID), Color_Channel_Structure}]); - disp('Error in tform_400, File will be saved, but the most crude matching was not made properly, double check the original image and process image'); - img_warped_no_scale = img_act_no_scale; -% end - end - - try - tform_noScale_200 = fitgeotrans(matchedPoints_act_200_tuned,matchedPoints_ref_200_tuned,'pwl'); - [ img_warped_no_scale,~] = imwarp_custom((img_warped_no_scale), tform_noScale_200,... - 'cubic','OutputView',outputView_no_scale,'WarpRescale',resc_warp); - if ~isempty(cell_warped_x) - distfun = @(p) norm(transformPointsInverse(tform_noScale_200, p) - double([cell_warped_x, cell_warped_y])); - options = optimset('Display','none','MaxFunEvals', numel(cell_warped_x)*500, 'MaxIter', numel(cell_warped_x)*500); - [cell_warpedXY, err] = fminsearch(distfun, double([cell_warped_x, cell_warped_y]), options); - cell_warped_x = cell_warpedXY(:,1); - cell_warped_y = cell_warpedXY(:,2); - end - catch -% disp([img_name{img_idx(img_ID), Color_Channel_Structure}, ': PWL fails in 200um scale - Trying LWM']); -% try -% tform_noScale_200 = fitgeotrans(matchedPoints_act_200_tuned,matchedPoints_ref_200_tuned,'lwm', 12); -% [ img_warped_no_scale,~] = imwarp_custom((img_warped_no_scale), tform_noScale_200,... -% 'cubic','OutputView',outputView_no_scale,'WarpRescale',resc_warp); -% if ~isempty(cell_warped_x) -% distfun = @(p) norm(transformPointsInverse(tform_noScale_200, p) - double([cell_warped_x, cell_warped_y])); -% options = optimset('Display','none','MaxFunEvals', numel(cell_warped_x)*500, 'MaxIter', numel(cell_warped_x)*500); -% [cell_warpedXY, err] = fminsearch(distfun, double([cell_warped_x, cell_warped_y]), options); -% cell_warped_x = cell_warpedXY(:,1); -% cell_warped_y = cell_warpedXY(:,2); -% end -% catch - disp(['PWL fails during analyzing : ', img_name{img_idx(img_ID), Color_Channel_Structure}]); - disp('File will be saved, but the 200um matching was not made properly, double check the original image and process image'); -% end - end - - try - tform_noScale_100 = fitgeotrans(matchedPoints_act_100_tuned,matchedPoints_ref_100_tuned,'pwl'); - [ img_warped_no_scale,~] = imwarp_custom(img_warped_no_scale, tform_noScale_100,... - 'cubic','OutputView',outputView_no_scale,'WarpRescale',resc_warp); - - if ~isempty(cell_warped_x) - distfun = @(p) norm(transformPointsInverse(tform_noScale_100, p) - double([cell_warped_x, cell_warped_y])); - options = optimset('Display','none','MaxFunEvals', numel(cell_warped_x)*500, 'MaxIter', numel(cell_warped_x)*500); - [cell_warpedXY, err] = fminsearch(distfun, double([cell_warped_x, cell_warped_y]), options); - cell_warped_x = cell_warpedXY(:,1); - cell_warped_y = cell_warpedXY(:,2); - end - catch -% disp([img_name{img_idx(img_ID), Color_Channel_Structure}, ': PWL fails in 100um scale - Trying LWM']); -% try -% tform_noScale_100 = fitgeotrans(matchedPoints_act_100_tuned,matchedPoints_ref_100_tuned,'lwm', 12); -% [ img_warped_no_scale,~] = imwarp_custom(img_warped_no_scale, tform_noScale_100,... -% 'cubic','OutputView',outputView_no_scale,'WarpRescale',resc_warp); -% -% if ~isempty(cell_warped_x) -% distfun = @(p) norm(transformPointsInverse(tform_noScale_100, p) - double([cell_warped_x, cell_warped_y])); -% options = optimset('Display','none','MaxFunEvals', numel(cell_warped_x)*500, 'MaxIter', numel(cell_warped_x)*500); -% [cell_warpedXY, err] = fminsearch(distfun, double([cell_warped_x, cell_warped_y]), options); -% cell_warped_x = cell_warpedXY(:,1); -% cell_warped_y = cell_warpedXY(:,2); -% end -% catch - disp(['PWL fails during analyzing : ', img_name{img_idx(img_ID), Color_Channel_Structure}]); - disp('File will be saved, but please check the image afterwards (100um warping was not made)'); -% end - end - - try - tform_noScale_50 = fitgeotrans(matchedPoints_act_50_tuned,matchedPoints_ref_50_tuned,'pwl'); - [ img_warped_no_scale,~] = imwarp_custom(img_warped_no_scale, tform_noScale_50,... - 'cubic','OutputView',outputView_no_scale,'WarpRescale',resc_warp); - - if ~isempty(cell_warped_x) - distfun = @(p) norm(transformPointsInverse(tform_noScale_50, p) - double([cell_warped_x, cell_warped_y])); - options = optimset('Display','none','MaxFunEvals', numel(cell_warped_x)*500, 'MaxIter', numel(cell_warped_x)*500); - [cell_warpedXY, err] = fminsearch(distfun, double([cell_warped_x, cell_warped_y]), options); - cell_warped_x = cell_warpedXY(:,1); - cell_warped_y = cell_warpedXY(:,2); - end - catch -% disp([img_name{img_idx(img_ID), Color_Channel_Structure}, ': PWL fails in 50um scale - Trying LWM']); -% try -% tform_noScale_50 = fitgeotrans(matchedPoints_act_50_tuned,matchedPoints_ref_50_tuned,'lwm', 12); -% [ img_warped_no_scale,~] = imwarp_custom(img_warped_no_scale, tform_noScale_50,... -% 'cubic','OutputView',outputView_no_scale,'WarpRescale',resc_warp); -% -% if ~isempty(cell_warped_x) -% distfun = @(p) norm(transformPointsInverse(tform_noScale_50, p) - double([cell_warped_x, cell_warped_y])); -% options = optimset('Display','none','MaxFunEvals', numel(cell_warped_x)*500, 'MaxIter', numel(cell_warped_x)*500); -% [cell_warpedXY, err] = fminsearch(distfun, double([cell_warped_x, cell_warped_y]), options); -% cell_warped_x = cell_warpedXY(:,1); -% cell_warped_y = cell_warpedXY(:,2); -% end -% catch - disp(['PWL fails during analyzing : ', img_name{img_idx(img_ID), Color_Channel_Structure}]); - disp('File will be saved, but please check the image afterwards (50um warping was not made)'); -% end - end - - %%% Match size of images - whos_bigger=size(img_ANO)-size(img_warped_no_scale); - if whos_bigger(1)>=0; x_dim=size(img_ANO,1); else; x_dim=size(img_warped_no_scale,1); end - if whos_bigger(2)>=0; y_dim=size(img_ANO,2); else; y_dim=size(img_warped_no_scale,2); end - - img_warped_sizeMatch=[x_dim,y_dim]-size(img_warped_no_scale); img_warped_sizeMatch(img_warped_sizeMatch<0)=0; - img_ANO_sizeMatch=[x_dim,y_dim]-size(img_ANO); img_ANO_sizeMatch(img_ANO_sizeMatch<0)=0; - - img_ANO=padarray(img_ANO,img_ANO_sizeMatch,'post'); - img_warped_no_scale=padarray(img_warped_no_scale,img_warped_sizeMatch,'post'); - - img_ch_thumbnail=imresize(img_warped_no_scale,xy_pix/ref_atlas_vox_res); - [injection_site_pts] = injection_volume_pix( img_ch_thumbnail ); - img_essence(img_ID).injection.Color(color_ch_ID).valid_pts=injection_site_pts; - img_essence(img_ID).injection.Color(color_ch_ID).image_thumbnail=img_ch_thumbnail; - img_essence(img_ID).Color_image_downscaled(color_ch_ID).img=img_warped_no_scale; - - %%%% for process monitoring %%%% - - monitoring_img=img_warped_no_scale; - monitoring_img_raw=monitoring_img; - - ROI_map=cat(3,false(size(img_ANO)),logical(img_ANO),logical(img_ANO)); - - monitoring_img=imadjust(monitoring_img,stretchlim(monitoring_img,0),[0 1]); - monitoring_img=cat(3,monitoring_img,monitoring_img,monitoring_img); - monitoring_img(ROI_map)=255; - - %%% detect cells across the whole slice image - cell_detected_all = round([cell_warped_x, cell_warped_y]); - - if ~isempty(cell_detected_all) - out_bnd_alpha = ref_boundarypad_0809_step5( img_warped_no_scale, xy_pix/xy_pix_resc_factor ); - out_bnd=inShape(out_bnd_alpha,cell_detected_all(:,2),cell_detected_all(:,1)); - cell_detected_all=cell_detected_all(~out_bnd,:); - end - - if ~isempty(cell_detected_all) - cell_detected_all_pos_ind=sub2ind(size(img_warped_no_scale),... - cell_detected_all(:,2),cell_detected_all(:,1)); - else - cell_detected_all_pos_ind=[]; - end - - img_essence(img_ID).Color_Cells(color_ch_ID).cell_locations=cell_detected_all; - - image_analyzed_ROI=figure; - imshow(monitoring_img,[]); hold on - title(strcat({'Image Name : '}, img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},... - {' - Ch No. '},num2str(color_ch_ID)),'Interpreter', 'none'); - if ~isempty(cell_detected_all) - scatter(cell_detected_all(:,1),cell_detected_all(:,2),9,'r','filled') - end + %%% Recalculate transform parameters for warping CELL COORDINATES %%% + transform_params_cellPos_downscaled = ... + invertElastixTransform(elastixDir_currImg); + %%%NOTE: I know the line above seems redundant and costs you + %%%extra few minutes for recalculation, but transforming the + %%%cell coordinates using "transform_params_ori_scale" (used for + %%%image warping) kept giving me weird results. It also doesn't + %%%work the other way round (i.e.img warping with + %%%"transform_params_cellPos_downscaled") + + img_essence(img_ID).transform_params_img_ori_scale= ... + transform_params_img_ori_scale; + img_essence(img_ID).transform_params_cellPos_downscaled= ... + transform_params_cellPos_downscaled; + + + %%% Load/Warp Images and Cell Pos in Original Resolution %%% + for color_ch_ID=1:length(Color_Channel_Interest) + cells_in_orig_img = ... + cell_detection_rs(img_ID).Color_Cells(color_ch_ID).cell_locations; + + img_Color=(imread(img_name{img_idx(img_ID),... + Color_Channel_Interest(color_ch_ID)})); + try; img_Color=rgb2gray(img_Color); end + img_Color=img_Color(img_info(img_idx(img_ID)).slice_window(1):... + img_info(img_idx(img_ID)).slice_window(2),... + img_info(img_idx(img_ID)).slice_window(3):... + img_info(img_idx(img_ID)).slice_window(4)); + img_Color=img_Color.*img_act_pad; + img_Color=padarray(img_Color,round([3000 3000]/(xy_pix))); + + + [img_warped_no_scale,~]=transformix(img_Color,... + transform_params_img_ori_scale); + img_warped_no_scale=uint8(img_warped_no_scale); + + [cell_detected_all,~]=transformix(cells_in_orig_img/mov_scale_factor, ... + transform_params_cellPos_downscaled); + cell_detected_all=round(cell_detected_all.OutputPoint*... + mov_scale_factor); %%% Looks a bit weird, but works.... + + %%% Match size of annotation and warped images + whos_bigger=size(img_ANO)-size(img_warped_no_scale); + if whos_bigger(1)>=0; x_dim=size(img_ANO,1); + else; x_dim=size(img_warped_no_scale,1); end + if whos_bigger(2)>=0; y_dim=size(img_ANO,2); + else; y_dim=size(img_warped_no_scale,2); end + + img_warped_sizeMatch=[x_dim,y_dim]-size(img_warped_no_scale); + img_warped_sizeMatch(img_warped_sizeMatch<0)=0; + img_ANO_sizeMatch=[x_dim,y_dim]-size(img_ANO); + img_ANO_sizeMatch(img_ANO_sizeMatch<0)=0; + img_ANO=padarray(img_ANO,img_ANO_sizeMatch,'post'); + img_warped_no_scale=padarray(img_warped_no_scale,... + img_warped_sizeMatch,'post'); + + %%% nope nope nope don't use the below... + % [injection_site_pts] = injection_volume_pix( img_ch_thumbnail ); + % img_essence(img_ID).injection.Color(color_ch_ID).valid_pts=injection_site_pts; + % img_essence(img_ID).injection.Color(color_ch_ID).image_thumbnail=img_ch_thumbnail; + % img_essence(img_ID).Color_image_downscaled(color_ch_ID).img=img_warped_no_scale; + + + %%%% for process monitoring (IMG) %%%% + monitoring_img=img_warped_no_scale; + monitoring_img_raw=monitoring_img; + + ROI_map=cat(3,... + false(size(img_ANO)),logical(img_ANO),logical(img_ANO)); + + monitoring_img=imadjust(monitoring_img,... + stretchlim(monitoring_img,0),[0 1]); + monitoring_img=cat(3,... + monitoring_img,monitoring_img,monitoring_img); + monitoring_img(ROI_map)=255; + + %%%% for process monitoring (CELL POS) %%%% + if ~isempty(cell_detected_all) + out_bnd_alpha = ref_boundarypad_0809_step5(... + img_warped_no_scale, xy_pix ); + out_bnd=inShape(out_bnd_alpha,... + cell_detected_all(:,2),cell_detected_all(:,1)); + cell_detected_all=cell_detected_all(~out_bnd,:); + end + + if ~isempty(cell_detected_all) + cell_detected_all_pos_ind=sub2ind(size(img_warped_no_scale),... + cell_detected_all(:,2),cell_detected_all(:,1)); + else + cell_detected_all_pos_ind=[]; + end + + img_essence(img_ID).Color_Cells(color_ch_ID).cell_locations=... + cell_detected_all; + + %%%% for process monitoring (All TOGETHER) %%%% + image_analyzed_ROI=figure; %%% with ROI boundaries overlaid + imshow(monitoring_img,[]); hold on + title(strcat({'Image Name : '}, ... + img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},... + {' - Ch No. '},num2str(color_ch_ID)),'Interpreter', 'none'); + if ~isempty(cell_detected_all) + scatter(cell_detected_all(:,1),... + cell_detected_all(:,2),9,'r','filled') + end + + image_analyzed=figure; %%% w/o ROI boundaries overlaid + imshow(monitoring_img_raw,[]); hold on + title(strcat({'Image Name : '}, ... + img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},... + {' - Ch No. '},num2str(color_ch_ID)),'Interpreter', 'none'); + if ~isempty(cell_detected_all) + scatter(cell_detected_all(:,1),... + cell_detected_all(:,2),9,'r','filled') + end + + try + save_name=strcat('/Image_Analysed/',... + img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},'.fig'); + saveas( image_analyzed_ROI ,[pwd save_name]); + close(image_analyzed_ROI) - image_analyzed=figure; - imshow(monitoring_img_raw,[]); hold on - title(strcat({'Image Name : '}, img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},... - {' - Ch No. '},num2str(color_ch_ID)),'Interpreter', 'none'); - if ~isempty(cell_detected_all) - scatter(cell_detected_all(:,1),cell_detected_all(:,2),9,'r','filled') - end - - try - save_name=strcat('/Image_Analysed/',img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},'.fig'); - saveas( image_analyzed_ROI ,[pwd save_name]); - close(image_analyzed_ROI) - - save_name=strcat('/Image_Analysed_ROI_absent/',img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},'.fig'); - saveas( image_analyzed,[pwd save_name]); - close(image_analyzed) - catch - disp('If this error message occurs without any error message above, something is wrong during saving stage, not from image processing stage'); - end + save_name=strcat('/Image_Analysed_ROI_absent/',... + img_name{img_idx(img_ID), Color_Channel_Interest(color_ch_ID)},'.fig'); + saveas( image_analyzed,[pwd save_name]); + close(image_analyzed) end end + end parfor_progress; end diff --git a/STEP_7_Reconstruction.m b/STEP_7_Reconstruction.m index 9f4fb40..53cb7bc 100644 --- a/STEP_7_Reconstruction.m +++ b/STEP_7_Reconstruction.m @@ -6,7 +6,7 @@ load Step5_ANO_info.mat; load ROI_region_pixels.mat; STEP_0_Parameters; -xy_pix=xy_pix/xy_pix_resc_factor; + %% Load Atlases [VOL, metaVOL] = nrrdread('ara_nissl_25_2017.nrrd');