ht_patch/apply_patch.pro at master · yyu1/ht_patch · GitHub

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
;Procedure for applying patching

;win_size is number of pixels in each direction from the center pixel: so win_size = 1 is 3x3

Pro apply_patch, in_ht_file, in_hv_file, out_ht_file, xdim, ydim, win_size

	COMMON setting, thresh   ;common block for values of settings

	thresh = 0.04  ; Use 0.04 for HV threshold

	max_ht = 65

	in_ht_image = fltarr(xdim,win_size*2+1)
	in_hv_image = fltarr(xdim,win_size*2+1)

	openr, in_ht_lun, in_ht_file, /get_lun
	openr, in_hv_lun, in_hv_file, /get_lun
	openw, out_lun, out_ht_file, /get_lun

	out_line = fltarr(xdim)
	hv_line = fltarr(xdim)

	ht_win = fltarr(win_size*2+1,win_size*2+1)
	hv_win = fltarr(win_size*2+1,win_size*2+1)

	print, 'Patching...'
	;copy first win_size lines, we patch but the pixel may not be the center of window
	readu, in_ht_lun, in_ht_image
	readu, in_hv_lun, in_hv_image

	for j=0ULL, win_size-1 do begin
		out_line[*] = in_ht_image[*,j]
		hv_line[*] = in_hv_image[*,j]

		index = where((out_line eq 0) and (hv_line ge thresh), count)

		if(count gt 0) then begin
			for i=0, count-1 do begin
				i_ind = index[i]
				i_ind_min = ((i_ind-win_size > 0) < (xdim-win_size*2-1))
				i_ind_max = ((i_ind+win_size > win_size*2) < (xdim-1))
				ht_win[*] = in_ht_image[i_ind_min:i_ind_max,*]
				hv_win[*] = in_hv_image[i_ind_min:i_ind_max,*]

				out_line[i_ind] = patch(ht_win,hv_win,hv_win[win_size,win_size])
			endfor
		endif
		writeu, out_lun, out_line < max_ht
		;writeu, out_lun, out_line
	endfor

	for j=ulong(win_size),ydim-win_size-1 do begin
		if (j mod 10000 eq 0) then print, j
		;set file pointer
		point_lun, in_ht_lun, (j-win_size)*xdim*4ULL
		point_lun, in_hv_lun, (j-win_size)*xdim*4ULL

		readu, in_ht_lun, in_ht_image
		readu, in_hv_lun, in_hv_image

		out_line[*] = in_ht_image[*,win_size]
		hv_line[*] = in_hv_image[*,win_size]

		index = where((out_line eq 0) and (hv_line ge thresh), count)

		if(count gt 0) then begin
			for i=0, count-1 do begin
				i_ind = index[i]
				i_ind_min = ((i_ind-win_size > 0) < (xdim-win_size*2-1))
				i_ind_max = ((i_ind+win_size > win_size*2) < (xdim-1))
				ht_win[*] = in_ht_image[i_ind_min:i_ind_max,*]
				hv_win[*] = in_hv_image[i_ind_min:i_ind_max,*]

				out_line[i_ind] = patch(ht_win,hv_win,hv_win[win_size,win_size])
			endfor
		endif

		writeu, out_lun, out_line < max_ht
		;writeu, out_lun, out_line
	endfor

	;patch the last win_size lines, since we read the block in the last step of previous loop, we don't need to read input anymore
	for j=0ULL, win_size-1 do begin
		out_line[*] = in_ht_image[*,win_size+j+1]
		hv_line[*] = in_hv_image[*,win_size+j+1]

		index = where((out_line eq 0) and (hv_line ge thresh), count)

		if(count gt 0) then begin
			for i=0, count-1 do begin
				i_ind = index[i]
				i_ind_min = ((i_ind-win_size > 0) < (xdim-win_size*2-1))
				i_ind_max = ((i_ind+win_size > win_size*2) < (xdim-1))
				ht_win[*] = in_ht_image[i_ind_min:i_ind_max,*]
				hv_win[*] = in_hv_image[i_ind_min:i_ind_max,*]

				out_line[i_ind] = patch(ht_win,hv_win,hv_win[win_size,win_size])
			endfor
		endif

		writeu, out_lun, out_line < max_ht
		;writeu, out_lun, out_line

	endfor

	free_lun, in_ht_lun, in_hv_lun, out_lun

End