-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathnew_pixel_sortv9.m
381 lines (324 loc) · 14 KB
/
new_pixel_sortv9.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
format long
clear all
close all
%folder = input('Assuming you are Carly on Colossus. Is this correct? (Y/N) ','s');
%if (folder == 'n') || (folder == 'N')
% folder = input('Please specify correct folder name/path: ','s');
%else
folder = '/home/derek/K2_PHunt/Data';
folder = '/media/derek/TOSHIBA/K2_data/K2_456/test';
folder = '/media/derek/TOSHIBA/K2_data/Dhital/test';
%folder = '/media/derek/TOSHIBA/temp';
%folder = '/media/derek/TOSHIBA/temp/temp2';
%%folder = '/media/derek/TOSHIBA/K2_data/PH';
%folder = '/media/derek/TOSHIBA/Clayton';
folder = '/media/derek/TOSHIBA/K2_data/Stello/C8/raw';
folder = '/media/derek/TOSHIBA/K2_data/Stello/C5';
folder = '/media/derek/TOSHIBA/K2_data/PH/2017';
folder = '/media/derek/TOSHIBA/K2_data/WideBinaries';
%folder = '/media/derek/TOSHIBA/K2_data/Trappist';
folder = '/media/derek/TOSHIBA/K2_data/PH/2017/C10';
%end
% uncommented for loop allows multiple files to be run at once:
% file = input('full file name: ','s');
% file = strcat(folder,'/',file); %
kplrfiles = dir(strcat(folder,'/ktwo*'));
for fileNum = 1:(length(kplrfiles))
%kplrfiles = dir(folder)
%for fileNum = 3:(length(kplrfiles))
file = kplrfiles(fileNum).name;
file = strcat(folder,'/',file);
% disp('getting data')
clear time flux bkg xcen ycen
clear mask
[time,apdim,n,dim,series,data,cnum] = get_k2_data(file);
series(isnan(series))= 0;
% meanseries = mean(series(:))
% disp('calling mean_image_gen')
sat_flag = 0;
[mean_image] = mean_image_gen(series,dim,n,file);
tmax = max(mean_image(:));
% if max(mean_image(:))>1.8e5
% sat_flag = 1;
% end
%Need to get smarter about how we handle saturated frames, to deal
%with the case where there's a saturated and a nonsaturated star in
%the same frame. I think the code handles it OK, the problem is
%identifying that a saturated or near-saturated cluster is really only ONE target and not
%multiples...
% disp('back from mean_image_gen')
filename = strrep(file,'_lpd-targ.fits','');
filename = strrep(filename,'ktwo','outputs/pipeout_ktwo');
filename = strcat(filename,'_','mean_image');
dlmwrite(filename, mean_image);
%Characterize and remove background
bkg = fit_background2(apdim,series);
% disp('Background removed')
for count = 1:apdim
series(:,:,count) = series(:,:,count) - bkg(count);
temp = series(:,:,count);
if min(temp(:))<0
series(:,:,count) = series(:,:,count)-min(temp(:));
end
end
%sat_flag = 1;
if sat_flag>0
%First start to handle saturated frames
test_image = zeros(dim(1),dim(2));
test_image(mean_image>3*mean(bkg))=1;
cc = bwconncomp(test_image);
numPixels = cellfun(@numel,cc.PixelIdxList);
[biggest,idx] = max(numPixels);
mask{1} = cc.PixelIdxList{idx};
num_mask = 1;
num = num_mask;
for count = 1:apdim
temp = series(:,:,count);
temp = temp(:);
test_image(~cc.PixelIdxList{idx}) = 0;
temp = temp.*test_image(:);
flux(1,count) = sum(temp(:));
end
%adjust mask-finding routine to be water-level type....but start at
%level where there's only one? or just a fractin of the peak value?
% and LOWER water level until there are >1 or until 3-bkg level
% is reached
%%OK probably should update variable names above
%%next extend to non-saturated case
else
% temp_image = mean_image;
temp_image = zeros(dim(1),dim(2));
cutoff = max(5*mean(bkg),2e-3*max(mean_image(:)));
temp_image(mean_image>cutoff)=1;
temp_image = mean_image;
cutoff = min(3*mean(bkg),max(mean_image(:))-1);
temp_image(mean_image<cutoff)=0; %undeleted this one
if max(mean_image(:))>1.8e5
temp_image = imgaussfilt(temp_image,1);
else
temp_image = imgaussfilt(temp_image);
end
bw = imregionalmax(temp_image);
[out,num] = bwlabel(bw);
%now we've determined the number num of maxima in the image
%remove the maxes that are too faint...
% stop
%%still need to fix what happens when num = 1, and also do we
%%still need the saturated kluge?
if num>1
tmp_num = num;
for i=1:num
tmp_out = out;
tmp_out(tmp_out~=i) = 0;
tmp_out = tmp_out/i;
max_test = mean_image.*tmp_out;
if max(max_test(:)) < 50 && tmp_num > 1 %is this how we want to set num?
tmp_num = tmp_num-1;
out(out==i) = 0;
end
end
num = tmp_num;
end
%% at this point the number of targets in the frame, num, has been determined
id_vals = unique(out);
if id_vals(1) == 0
id_vals(1) = [];
end
clear peak_val
for i = 1:numel(id_vals)
peak_val(i) = mean_image(out == id_vals(i)); %%these are the values at the peaks
end
%% at this point the values of the peaks for each target have been determined
%% adding a section to handle situation where there is only one mask
if num == 1
%divide range between peak and mean background level into
%100 pieces
min_val = ceil(mean(bkg(:)));
%min_val = 1000;
max_val = floor(peak_val(1))-1;
range_val = 0.01*(max_val-min_val);
%for each value need to calculate the relevant statistic...
for calc_val = 1:100
t_image = mean_image;
t_image(t_image<(min_val+calc_val*range_val)) = 0; %zero all values below the reference value
cct = bwconncomp(t_image,4);
mask{1} = cct.PixelIdxList{1};
clear tflux
for count = 1:apdim
temp = series(:,:,count);
temp = temp(mask{1});
tflux(count) = sum(temp(:));
end
%tstat(calc_val) = sum(abs(diff(diff(diff(tflux)))));
tstat(calc_val) = sum(abs(diff(tflux)));
tstat(calc_val) = tstat(calc_val)/mean(tflux);
%plot(tflux/mean(tflux))
%calc_val, std(tflux)/mean(tflux)
%hold on
end
[aa bb] = min(tstat);
t_image = mean_image;
t_image(t_image<(min_val+bb*range_val)) = 0;
%bb=1; %this approach doesn't work for saturated targets!
%%look at brightest 9 pixel group mean
sort_pix = sort(mean_image(:),'descend');
sort_pix = mean(sort_pix(1:9));
if sort_pix > 1.7e5
bb = 1;
end
%t_image(t_image<(min_val+bb*range_val)) = 0;
%t_image(t_image<450) = 0;
cct = bwconncomp(t_image);
%mask{1} = cct.PixelIdxList{2};
mask{1} = cct.PixelIdxList{1};
%numel(mask{1})
else
temp = find(out);
fin_ref_value = zeros(1,num);
for ii=1:num
ploc = temp(ii); % location of this peak
parr = setdiff(temp,ploc); %all peaks not this one
ref_value = ceil(peak_val(ii))-1; %start at value at this peak
%fin_ref_value(ii) = 1;
%check if current peak is inside
for j = ref_value:-1:0 %Good for step 1, but will have to replace with binary search tree
t_image = mean_image;
t_image(t_image<j) = 0; %zero all values below the reference value
cct = bwconncomp(t_image);
for i=1:cct.NumObjects
if ismember(ploc,cct.PixelIdxList{i}) & any(ismember(parr,cct.PixelIdxList{i}))
fin_ref_value(ii) = j+1; %overlapping masks...increase ref value and solution is found
%cct.PixelIdxList{ii} %Need to figure out how to label
%the elements in the mask!
end
end
if fin_ref_value(ii) > 0 %if this has been set, then we are done!
t_image = mean_image;
t_image(t_image<fin_ref_value(ii)) = 0;
cct = bwconncomp(t_image,4);
for k = 1:cct.NumObjects
% disp(ploc)
if ismember(ploc,cct.PixelIdxList{k})
% cct.PixelIdxList{1};
% cct.PixelIdxList{2};
mask{ii} = cct.PixelIdxList{k};
end
end
break
end
end
end
num = numel(mask);
%%I think in here is where we need to check mask sizes for the multi-mask
%%case, to see if shrinking the mask improves the SNR of the dewhitened
%%time series. Should be able to steal code from single-mask case!
%for ii=1:num
% %divide range between peak and mean background level into
% %100 pieces
% min_val = ceil(fin_ref_value(ii));
% %min_val = 1000;
% max_val = floor(peak_val(ii))-1;
% range_val = 0.01*(max_val-min_val);
% %for each value need to calculate the relevant statistic...
% for calc_val = 1:100
% t_image = mean_image;
% t_image(t_image<(min_val+calc_val*range_val)) = 0; %zero all values below the reference value
% cct = bwconncomp(t_image);
% mask{ii} = cct.PixelIdxList{ii};
% for count = 1:apdim
% temp = series(:,:,count);
% temp = temp(mask{ii});
% tflux(count) = sum(temp(:));
% end
% tstat(calc_val) = sum(abs(diff(diff(diff(tflux)))));
% tstat(calc_val) = tstat(calc_val)/mean(tflux);
% end
%
% [aa bb] = min(tstat);
% t_image = mean_image;
% t_image(t_image<(min_val+bb*range_val)) = 0;
% cct = bwconncomp(t_image);
% mask{ii} = cct.PixelIdxList{ii};
%end
end
%% at this point the masks have been determined
end
% end
%Calculate fluxes and centroids
for i=1:num
%prep for centroid calculations
%first convert to x-y positions
[xx yy] = ind2sub([dim(1), dim(2)],mask{i});
for count = 1:apdim
temp = series(:,:,count);
temp = temp(mask{i});
flux(i,count) = sum(temp(:));
%calculate
xflux = 0.0;
yflux = 0.0;
for j = 1:numel(mask{i})
xflux = xflux + xx(j)*series(xx(j),yy(j),count);
yflux = yflux + yy(j)*series(xx(j),yy(j),count);
end
xcen(i,count) = xflux/flux(i,count);
ycen(i,count) = yflux/flux(i,count);
end
xc(i) = mean(xcen(i,:));
yc(i) = mean(ycen(i,:));
end
%%temporarily try centroids which are simply the weighted sum of the entire background-subtracted frame
for i=1:num
%prep for centroid calculations
%first convert to x-y positions
[xx yy] = ind2sub([dim(1), dim(2)],[1:1:dim(1)*dim(2)]);
for count = 1:apdim
temp = series(:,:,count);
% temp = temp(mask{i});
tflux(i,count) = sum(temp(:));
%calculate
xflux = 0.0;
yflux = 0.0;
for j = 1:(dim(1)*dim(2))
xflux = xflux + xx(j)*series(xx(j),yy(j),count);
yflux = yflux + yy(j)*series(xx(j),yy(j),count);
end
xcen(i,count) = xflux/tflux(i,count);
ycen(i,count) = yflux/tflux(i,count);
end
xcen(i,:) = xcen(i,:)-mean(xcen(i,:))+xc(i);
ycen(i,:) = ycen(i,:)-mean(ycen(i,:))+yc(i);
end
%reorder based on mean flux!
% reorganize mask order so 1 is brightest
clear nmask
clear nmean
nmask = mask;
nflux = flux;
nxcen = xcen;
nycen = ycen;
if num > 1
for i=1:num
nmean(i) = mean(flux(i,:));
end
[ai bi] = sort(nmean,'descend');
for i=1:num
mask{i} = nmask{bi(i)};
flux(i,:) = nflux(bi(i),:);
xcen(i,:) = nxcen(bi(i),:);
ycen(i,:) = nycen(bi(i),:);
end
end
%put raw lightcurves out to file
out_to_file2(time,flux,bkg,xcen,ycen,file,num);
% export image of mask
mask_image = zeros(dim(1),dim(2));
for i=1:num
mask_image(mask{i}) = i;
end
hmask = imagesc(mask_image);
colorbar;
filename = strrep(file,'_lpd-targ.fits','');
filename = strrep(filename,'ktwo','outputs/pipeout_ktwo');
filename = strcat(filename,'_','mask.jpg');
saveas(hmask,filename);
end