-
Notifications
You must be signed in to change notification settings - Fork 0
/
unit_stack_routines.pas
1829 lines (1576 loc) · 95.7 KB
/
unit_stack_routines.pas
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
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
unit unit_stack_routines;
{Copyright (C) 2017, 2021 by Han Kleijn, www.hnsky.org
email: han.k.. at...hnsky.org
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at https://mozilla.org/MPL/2.0/. }
{$mode delphi}
interface
uses
Classes, SysUtils,forms, math, unit_stack, astap_main, unit_star_align;
procedure stack_LRGB(oversize:integer; var files_to_process : array of TfileToDo; out counter : integer );{stack LRGB mode}
procedure stack_average(oversize:integer; var files_to_process : array of TfileToDo; out counter : integer);{stack average}
procedure stack_mosaic(oversize:integer; var files_to_process : array of TfileToDo;max_dev_backgr: double; out counter : integer);{stack mosaic/tile mode}
procedure stack_sigmaclip(oversize:integer; var files_to_process : array of TfileToDo; out counter : integer); {stack using sigma clip average}
procedure calibration_and_alignment(oversize:integer; var files_to_process : array of TfileToDo; out counter : integer); {calibration_and_alignment only}
{$inline on} {!!! Set this off for debugging}
procedure calc_newx_newy(vector_based : boolean; fitsXfloat,fitsYfloat: double); inline; {apply either vector or astrometric correction}
procedure astrometric_to_vector; {convert astrometric solution to vector solution}
procedure initialise_var1;{set variables correct}
procedure initialise_var2;{set variables correct}
function test_bayer_matrix(img: image_array) :boolean; {test statistical if image has a bayer matrix. Execution time about 1ms for 3040x2016 image}
var
pedestal_s : double;{target background value}
var
SIN_dec0,COS_dec0,x_new_float,y_new_float,SIN_dec_ref,COS_dec_ref,
ap_0_1_ref,ap_0_2_ref,ap_0_3_ref,ap_1_0_ref,ap_1_1_ref, ap_1_2_ref,ap_2_0_ref,ap_2_1_ref,ap_3_0_ref, bp_0_1_ref,bp_0_2_ref,bp_0_3_ref,bp_1_0_ref,bp_1_1_ref,bp_1_2_ref,bp_2_0_ref,bp_2_1_ref,bp_3_0_ref : double;
gain_refxxx: string;
implementation
uses unit_astrometric_solving;
procedure calc_newx_newy(vector_based : boolean; fitsXfloat,fitsYfloat: double); inline; {apply either vector or astrometric correction}
var
u,u0,v,v0,dRa,dDec,delta,ra_new,dec_new,delta_ra,det,gamma,SIN_dec_new,COS_dec_new,SIN_delta_ra,COS_delta_ra,h: double;
Begin
if vector_based then {vector based correction}
begin
x_new_float:=solution_vectorX[0]*(fitsxfloat-1)+solution_vectorX[1]*(fitsYfloat-1)+solution_vectorX[2]; {correction x:=aX+bY+c x_new_float in image array range 0..head.width-1}
y_new_float:=solution_vectorY[0]*(fitsxfloat-1)+solution_vectorY[1]*(fitsYfloat-1)+solution_vectorY[2]; {correction y:=aX+bY+c}
end
else
begin {astrometric based correction}
{6. Conversion (x,y) -> (RA,DEC) for image to be added}
u0:=fitsXfloat-head.crpix1;
v0:=fitsYfloat-head.crpix2;
if a_order>=2 then {apply SIP correction up second order}
begin
u:=u0 + a_2_0*u0*u0 + a_0_2*v0*v0 + a_1_1*u0*v0; {SIP correction}
v:=v0 + b_2_0*u0*u0 + b_0_2*v0*v0 + b_1_1*u0*v0; {SIP correction}
end
else
begin
u:=u0;
v:=v0;
end;
dRa :=(head.cd1_1 * u +head.cd1_2 * v)*pi/180;
dDec:=(head.cd2_1 * u +head.cd2_2 * v)*pi/180;
delta:=COS_dec0 - dDec*SIN_dec0;
gamma:=sqrt(dRa*dRa+delta*delta);
RA_new:=head.ra0+arctan(Dra/delta);
DEC_new:=arctan((SIN_dec0+dDec*COS_dec0)/gamma);
{5. Conversion (RA,DEC) -> (x,y) of reference image}
sincos(dec_new,SIN_dec_new,COS_dec_new);{sincos is faster then separate sin and cos functions}
delta_ra:=RA_new-head_ref.ra0;
sincos(delta_ra,SIN_delta_ra,COS_delta_ra);
H := SIN_dec_new*sin_dec_ref + COS_dec_new*COS_dec_ref*COS_delta_ra;
dRA := (COS_dec_new*SIN_delta_ra / H)*180/pi;
dDEC:= ((SIN_dec_new*COS_dec_ref - COS_dec_new*SIN_dec_ref*COS_delta_ra ) / H)*180/pi;
det:=head_ref.CD2_2*head_ref.CD1_1 - head_ref.CD1_2*head_ref.CD2_1;
u0:= - (head_ref.CD1_2*dDEC - head_ref.CD2_2*dRA) / det;
v0:= + (head_ref.CD1_1*dDEC - head_ref.CD2_1*dRA) / det;
if ap_order>=2 then {apply SIP correction up to second order}
begin
x_new_float:=(head_ref.crpix1 + u0+ap_0_1*v0+ ap_0_2*v0*v0+ ap_0_3*v0*v0*v0 +ap_1_0*u0 + ap_1_1*u0*v0+ ap_1_2*u0*v0*v0+ ap_2_0*u0*u0 + ap_2_1*u0*u0*v0+ ap_3_0*u0*u0*u0)-1;{3th order SIP correction, fits count from 1, image from zero therefore subtract 1}
y_new_float:=(head_ref.crpix2 + v0+bp_0_1*v0+ bp_0_2*v0*v0+ bp_0_3*v0*v0*v0 +bp_1_0*u0 + bp_1_1*u0*v0+ bp_1_2*u0*v0*v0+ bp_2_0*u0*u0 + bp_2_1*u0*u0*v0+ bp_3_0*u0*u0*u0)-1;{3th order SIP correction}
end
else
begin
x_new_float:=(head_ref.crpix1 + u0)-1; {in image array range 0..width-1}
y_new_float:=(head_ref.crpix2 + v0)-1;
end;
end;{astrometric}
end;{calc_newx_newy}
procedure astrometric_to_vector;{convert astrometric solution to vector solution}
var
flipped,flipped_reference : boolean;
begin
a_order:=0; {SIP correction should be zero by definition}
calc_newx_newy(false,1,1) ;
solution_vectorX[2]:=x_new_float;
solution_vectorY[2]:=y_new_float;
calc_newx_newy(false,2, 1); {move one pixel in X}
solution_vectorX[0]:=+(x_new_float- solution_vectorX[2]);
solution_vectorX[1]:=-(y_new_float- solution_vectorY[2]);
calc_newx_newy(false,1, 2);{move one pixel in Y}
solution_vectorY[0]:=-(x_new_float- solution_vectorX[2]);
solution_vectorY[1]:=+(y_new_float- solution_vectorY[2]);
flipped:=(head.cd1_1/head.cd2_2)>0; {flipped image. Either flipped vertical or horizontal but not both. Flipped both horizontal and vertical is equal to 180 degrees rotation and is not seen as flipped}
flipped_reference:=(head_ref.CD1_1/head_ref.CD2_2)>0; {flipped reference image}
if flipped<>flipped_reference then {this can happen is user try to add images from a diffent camera/setup}
begin
solution_vectorX[1]:=-solution_vectorX[1];
solution_vectorY[0]:=-solution_vectorY[0];
end;
end;
procedure initialise_var1;{set variables correct}
begin
sincos(head_ref.dec0,SIN_dec_ref,COS_dec_ref);{do this in advance since it is for each pixel the same}
end;
procedure initialise_var2;{set variables correct}
begin
ap_0_1_ref:=ap_0_1;{store polynomial first fits }
ap_0_2_ref:=ap_0_2;
ap_0_3_ref:=ap_0_3;
ap_1_0_ref:=ap_1_0;
ap_1_1_ref:=ap_1_1;
ap_1_2_ref:=ap_1_2;
ap_2_0_ref:=ap_2_0;
ap_2_1_ref:=ap_2_1;
ap_3_0_ref:=ap_3_0;
bp_0_1_ref:=bp_0_1;
bp_0_2_ref:=bp_0_2;
bp_0_3_ref:=bp_0_3;
bp_1_0_ref:=bp_1_0;
bp_1_1_ref:=bp_1_1;
bp_2_1_ref:=bp_2_1;
bp_2_0_ref:=bp_2_0;
bp_2_1_ref:=bp_2_1;
bp_3_0_ref:=bp_3_0;
end;
procedure stack_LRGB(oversize:integer; var files_to_process : array of TfileToDo; out counter : integer );{stack LRGB mode}
var
fitsX,fitsY,c,width_max, height_max, x_new,y_new, binning,oversizeV : integer;
background_r, background_g, background_b, background_l ,
rgbsum,red_f,green_f,blue_f, value ,colr, colg,colb, red_add,green_add,blue_add,
rr_factor, rg_factor, rb_factor,
gr_factor, gg_factor, gb_factor,
br_factor, bg_factor, bb_factor,
saturated_level,hfd_min : double;
init, solution,use_star_alignment,use_manual_align,use_ephemeris_alignment,
use_astrometry_internal,vector_based :boolean;
warning : string;
begin
with stackmenu1 do
begin
{move often uses setting to booleans. Great speed improved if use in a loop and read many times}
use_star_alignment:=stackmenu1.use_star_alignment1.checked;
use_manual_align:=stackmenu1.use_manual_alignment1.checked;
use_ephemeris_alignment:=stackmenu1.use_ephemeris_alignment1.checked;
use_astrometry_internal:=use_astrometry_internal1.checked;
hfd_min:=max(0.8 {two pixels},strtofloat2(stackmenu1.min_star_size_stacking1.caption){hfd});{to ignore hot pixels which are too small}
counter:=0;
jd_sum:=0;{sum of Julian midpoints}
jd_stop:=0;{end observations in Julian day}
init:=false;
{LRGB method}
begin
memo2_message('Combining colours.');
rr_factor:=strtofloat2(rr1.text);
rg_factor:=strtofloat2(rg1.text);
rb_factor:=strtofloat2(rb1.text);
gr_factor:=strtofloat2(gr1.text);
gg_factor:=strtofloat2(gg1.text);
gb_factor:=strtofloat2(gb1.text);
br_factor:=strtofloat2(br1.text);
bg_factor:=strtofloat2(bg1.text);
bb_factor:=strtofloat2(bb1.text);
background_r:=0;
background_g:=0;
background_b:=0;
background_l:=0;
red_add:=strtofloat2(red_filter_add1.text);
green_add:=strtofloat2(green_filter_add1.text);
blue_add:=strtofloat2(blue_filter_add1.text);
for c:=0 to length(files_to_process)-1 do {should contain reference,r,g,b,rgb,l}
begin
if c=5 then {all colour files added, correct for the number of pixel values added at one pixel. This can also happen if one colour has an angle and two pixel fit in one!!}
begin {fix RGB stack}
memo2_message('Applying black spot filter on interim RGB image.');
black_spot_filter(img_average);
end;{c=5, all colour files added}
if length(files_to_process[c].name)>0 then
begin
try { Do some lengthy operation }
filename2:=files_to_process[c].name;
if c=0 then memo2_message('Loading reference image: "'+filename2+'".');
if c=1 then memo2_message('Adding red file: "'+filename2+'" to final image.');
if c=2 then memo2_message('Adding green file: "'+filename2+'" to final image.');
if c=3 then memo2_message('Adding blue file: "'+filename2+'" to final image.');
if c=4 then memo2_message('Adding RGB file: "'+filename2+'" to final image.');
if c=5 then memo2_message('Using luminance file: "'+filename2+'" for final image.');
{load image}
Application.ProcessMessages;
if ((esc_pressed) or (load_fits(filename2,true {light}, true,init=false{update memo1} ,0,head,img_loaded)=false)) then begin memo2_message('Error');exit;end;{update memo for case esc is pressed}
if init=false then
begin
head_ref:=head;{backup solution}
initialise_var1;{set variables correct, do this before apply dark}
initialise_var2;{set variables correct}
end;
saturated_level:=head.datamax_org*0.97;{130}
if c=1 then
begin
get_background(0,img_loaded,true,false, {var} background_r,star_level);{unknown, do not calculate noise_level}
background_r:=background_r-500;{pedestal 500}
cblack:=round(background_r);
counterR:=head.light_count ;counterRdark:=head.dark_count; counterRflat:=head.flat_count; counterRbias:=head.flatdark_count; exposureR:=round(head.exposure);temperatureR:=head.set_temperature;{for historical reasons}
end;
if c=2 then
begin
get_background(0,img_loaded,true,false, {var} background_g,star_level);{unknown, do not calculate noise_level}
background_g:=background_g-500; {pedestal +500}
cblack:=round(background_g);
counterG:=head.light_count;counterGdark:=head.dark_count; counterGflat:=head.flat_count; counterGbias:=head.flatdark_count; exposureG:=round(head.exposure);temperatureG:=head.set_temperature;
end;
if c=3 then
begin
get_background(0,img_loaded,true,false, {var} background_b,star_level);{unknown, do not calculate noise_level}
background_b:=background_b-500; {pedestal +500}
cblack:=round( background_b);
counterB:=head.light_count; counterBdark:=head.dark_count; counterBflat:=head.flat_count; counterBbias:=head.flatdark_count; exposureB:=round(head.exposure);temperatureB:=head.set_temperature;
end;
if c=4 then
begin
get_background(0,img_loaded,true,false, {var} background_r,star_level);{unknown, do not calculate noise_level}
background_r:=background_r-500; {pedestal +500}
background_g:=background_r-500;
background_b:=background_r-500;
cblack:=round( background_r);
counterRGB:=head.light_count; counterRGBdark:=head.dark_count; counterRGBflat:=head.flat_count; counterRGBbias:=head.flatdark_count; exposureRGB:=round(head.exposure);;temperatureRGB:=head.set_temperature;
end;
if c=5 then {Luminance}
begin
get_background(0,img_loaded,true,false, {var} background_L,star_level);{unknown, do not calculate noise_level}
background_L:=background_L-500;
cblack:=round( background_L);
counterL:=head.light_count; counterLdark:=head.dark_count; counterLflat:=head.flat_count; counterLbias:=head.flatdark_count; exposureL:=round(head.exposure);temperatureL:=head.set_temperature;
end;
if use_astrometry_internal then {internal solver, create new solutions for the R, G, B and L stacked images if required}
begin
memo2_message('Preparing astrometric solution for interim file: '+filename2);
if head.cd1_1=0 then solution:= create_internal_solution(img_loaded) else solution:=true;
if solution=false {load astrometry.net solution succesfull} then begin memo2_message('Abort, No astrometric solution for '+filename2); exit;end;{no solution found}
//if c=0 then cdelt2_lum:=head.cdelt2;
//if ((c=1) or (c=2) or (c=3) or (c=4)) then
//begin
// if head.cdelt2>cdelt2_lum*1.01 then memo2_message('█ █ █ █ █ █ Warning!! RGB images have a larger pixel size ('+floattostrF(head.cdelt2*3600,ffgeneral,3,2)+'arcsec) then the luminance image ('+floattostrF(cdelt2_lum*3600,ffgeneral,3,2)+'arcsec). If the stack result is poor, select in tab "stack method" the 3x3 mean or 5x5 mean filter for smoothing interim RGB.');
//end;
end
else
if init=false then {first image}
begin
if ((use_manual_align) or (use_ephemeris_alignment)) then
begin
referenceX:=strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_X]); {reference offset}
referenceY:=strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_Y]); {reference offset}
end
else
begin
binning:=report_binning(head.height);{select binning based on the height of the light}
bin_and_find_stars(img_loaded, binning,1 {cropping},hfd_min,true{update hist},starlist1,warning);{bin, measure background, find stars}
find_quads(starlist1,0, quad_smallest,quad_star_distances1);{find quads for reference image/database}
end;
end;
if init=false then {init}
begin
if oversize<0 then {shrink a lot, adapt in ratio}
begin
oversize:=max(oversize,-round((head.width-100)/2) );{minimum image width is 100}
oversizeV:=round(oversize*head.height/head.width);{vertical shrinkage in pixels}
height_max:=head.height+oversizeV*2;
end
else
begin
oversizeV:=oversize;
height_max:=head.height+oversize*2;
end;
width_max:=head.width+oversize*2;
setlength(img_average,3,width_max,height_max);{will be color}
for fitsY:=0 to height_max-1 do
for fitsX:=0 to width_max-1 do
begin
img_average[0,fitsX,fitsY]:=0; {clear img_average}
img_average[1,fitsX,fitsY]:=0; {clear img_average}
img_average[2,fitsX,fitsY]:=0; {clear img_average}
end;
end;{init, c=0}
solution:=true;{assume solution is found}
if use_astrometry_internal then sincos(head.dec0,SIN_dec0,COS_dec0) {do this in advance since it is for each pixel the same}
else
begin {align using star match}
if init=true then {second image}
begin
if ((use_manual_align) or (use_ephemeris_alignment)) then
begin {manual alignment}
solution_vectorX[2]:=referenceX-strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_X]); {calculate correction}
solution_vectorY[2]:=referenceY-strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_Y]);
memo2_message('Solution x:=x+'+floattostr6(solution_vectorX[2])+', y:=y+'+floattostr6(solution_vectorY[2]));
end
else
begin{internal alignment}
bin_and_find_stars(img_loaded, binning,1 {cropping},hfd_min,true{update hist},starlist2,warning);{bin, measure background, find stars}
find_quads(starlist2,0, quad_smallest,quad_star_distances2);{find star quads for new image}
if find_offset_and_rotation(3,strtofloat2(stackmenu1.quad_tolerance1.text)) then {find difference between ref image and new image}
memo2_message(inttostr(nr_references)+' of '+ inttostr(nr_references2)+' quads selected matching within '+stackmenu1.quad_tolerance1.text+' tolerance.'
+' Solution x:='+floattostr6(solution_vectorX[0])+'*x+ '+floattostr6(solution_vectorX[1])+'*y+ '+floattostr6(solution_vectorX[2])
+', y:='+floattostr6(solution_vectorY[0])+'*x+ '+floattostr6(solution_vectorY[1])+'*y+ '+floattostr6(solution_vectorY[2]) )
else
begin
memo2_message('Not enough quad matches <3 or inconsistent solution, skipping this image.');
files_to_process[c].name:=''; {remove file from list}
solution:=false;
ListView1.Items.item[files_to_process[c].listviewindex].SubitemImages[L_result]:=6;{mark 3th column with exclaimation}
ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_result]:='no solution';{no stack result}
end;
end;{internal alignment}
end
else
reset_solution_vectors(1);{no influence on the first image}
end;{using star match}
init:=true;{initialize for first image done}
if ((c<>0) and (solution)) then {do not add reference channel c=0, in most case luminance file.}
begin
inc(counter);{count number of colour files involved}
date_to_jd(head.date_obs,head.exposure);{convert head.date_obs string and head.exposure time to global variables jd_start (julian day start head.exposure) and jd_mid (julian day middle of the head.exposure)}
if jd_start>jd_stop then jd_stop:=jd_start;{find latest start time}
jd_sum:=jd_sum+jd_mid;{sum julian days of images at midpoint head.exposure}
vector_based:=((use_star_alignment) or (use_manual_align) or (use_ephemeris_alignment));
if ((vector_based=false) and (a_order=0)) then {no SIP from astronomy.net}
begin
astrometric_to_vector;{convert astrometric solution to vector solution}
vector_based:=true;
end;
for fitsY:=1 to head.height do {skip outside pixels if color}
for fitsX:=1 to head.width do
begin
calc_newx_newy(vector_based,fitsX,fitsY);{apply correction}
x_new:=round(x_new_float+oversize);y_new:=round(y_new_float+oversizeV);
if ((x_new>=0) and (x_new<=width_max-1) and (y_new>=0) and (y_new<=height_max-1)) then
begin
if c=1 {red} then
begin
value:=img_loaded[0,fitsX-1,fitsY-1];
if value>saturated_level then {saturation, mark all three colors as black spot (<=0) to maintain star colour}
begin
img_average[0,x_new,y_new]:=0;//saturation marker, process later as black spot
img_average[1,x_new,y_new]:=0;
img_average[2,x_new,y_new]:=0;
end
else
begin
value:=(value-background_r);{image loaded is already corrected with dark and flat. Normalize background to level 500}{NOTE: fits count from 1, image from zero}
if rr_factor>0.00001 then begin img_average[0,x_new,y_new]:=rr_factor*value;{execute only if greater then zero for speed} end;
if rg_factor>0.00001 then begin img_average[1,x_new,y_new]:=rg_factor*value; end;
if rb_factor>0.00001 then begin img_average[2,x_new,y_new]:=rb_factor*value; end;
end;
end;
if c=2 {green} then
begin
value:=img_loaded[0,fitsX-1,fitsY-1];
if value>saturated_level then {saturation, mark all three colors as black spot (<=0) to maintain star colour}
begin
img_average[0,x_new,y_new]:=0;//saturation marker, process later as black spot
img_average[1,x_new,y_new]:=0;
img_average[2,x_new,y_new]:=0;
end
else
begin
value:=(value-background_g);{image loaded is already corrected with dark and flat. Normalize background to level 500}{NOTE: fits count from 1, image from zero}
if gr_factor>0.00001 then begin img_average[0,x_new,y_new]:=gr_factor*value;{execute only if greater then zero for speed} end;
if gg_factor>0.00001 then begin img_average[1,x_new,y_new]:=gg_factor*value; end;
if gb_factor>0.00001 then begin img_average[2,x_new,y_new]:=gb_factor*value; end;
end;
end;
if c=3 {blue} then
begin
value:=img_loaded[0,fitsX-1,fitsY-1];
if value>saturated_level then {saturation, mark all three colors as black spot (<=0) to maintain star colour}
begin
img_average[0,x_new,y_new]:=0;//saturation marker, process later as black spot
img_average[1,x_new,y_new]:=0;
img_average[2,x_new,y_new]:=0;
end
else
begin
value:=(value-background_b);{image loaded is already corrected with dark and flat. Normalize background to level 500}{NOTE: fits count from 1, image from zero}
if br_factor>0.00001 then begin img_average[0,x_new,y_new]:=br_factor*value;{execute only if greater then zero for speed} end;
if bg_factor>0.00001 then begin img_average[1,x_new,y_new]:=bg_factor*value; end;
if bb_factor>0.00001 then begin img_average[2,x_new,y_new]:=bb_factor*value; end;
end;
end;
if c=4 {RGB image, head.naxis3=3} then
begin
begin img_average[0,x_new,y_new]:=img_loaded[0,fitsX-1,fitsY-1]-background_r; end;
begin img_average[1,x_new,y_new]:=img_loaded[1,fitsX-1,fitsY-1]-background_g; end;
begin img_average[2,x_new,y_new]:=img_loaded[2,fitsX-1,fitsY-1]-background_b; end;
end;
if c=5 {Luminance} then
begin
{rgb is already blurred}
{r:=l*(0.33+r)/(r+g+b)}
colr:=img_average[0,x_new,y_new] - 475 + red_add; {lowest_most_common is around 450 to 500}
colg:=img_average[1,x_new,y_new] - 475 + green_add;
colb:=img_average[2,x_new,y_new] - 475 + blue_add;
rgbsum:=colr+colg+colb;
if rgbsum<0.1 then begin rgbsum:=0.1; red_f:=rgbsum/3; green_f:=red_f; blue_f:=red_f;end
else
begin
red_f:=colr/rgbsum; if red_f<0 then red_f:=0; if red_f>1 then red_f:=1;
green_f:=colg/rgbsum; if green_f<0 then green_f:=0;if green_f>1 then green_f:=1;
blue_f:=colb/rgbsum; if blue_f<0 then blue_f:=0; if blue_f>1 then blue_f:=1;
end;
img_average[0,x_new,y_new]:=1000+(img_loaded[0,fitsX-1,fitsY-1] - background_l)*(red_f);
img_average[1,x_new,y_new]:=1000+(img_loaded[0,fitsX-1,fitsY-1] - background_l)*(green_f);
img_average[2,x_new,y_new]:=1000+(img_loaded[0,fitsX-1,fitsY-1] - background_l)*(blue_f);
end;
end;
end;
end;
progress_indicator(94+c,' LRGB');{show progress, 95..99}
except
beep;
end;{try}
end;
end;
if counter<>0 then
begin
head:=head_ref; {restore solution. Works only if no oversize is used}
head.naxis3:=3;{three colours}
head.naxis :=3;{three dimensions. Header will be updated in the save routine}
img_loaded:=img_average;
head.width:=width_max;
head.height:=height_max;
end;
end;{LRGB}
end;{with stackmenu1}
end;
function test_bayer_matrix(img: image_array) :boolean; {test statistical if image has a bayer matrix. Execution time about 1ms for 3040x2016 image}
var
fitsX,w,h,middleY,step_size : integer;
p11,p12,p21,p22 : array of double;
m11,m12,m21,m22,lowest,highest : double;
const
steps=100;
begin
// colors:=Length(img); {colors}
w:=Length(img[0]); {width}
h:=Length(img[0][0]); {height}
middleY:=h div 2;
step_size:=w div steps;
if odd(step_size) then step_size:=step_size-1;{make even so it ends up at the correct location of the 2x2 matrix}
SetLength(p11,steps);
SetLength(p12,steps);
SetLength(p21,steps);
SetLength(p22,steps);
for fitsX:=0 to steps-1 do {test one horizontal line and take 100 samples of the bayer matrix}
begin
p11[fitsX]:=img[0,step_size*fitsX,middleY];
p12[fitsX]:=img[0,step_size*fitsX+1,middleY];
p21[fitsX]:=img[0,step_size*fitsX,middleY+1];
p22[fitsX]:=img[0,step_size*fitsX+1,middleY+1];
end;
m11:=Smedian(p11,steps);
m12:=Smedian(p12,steps);
m21:=Smedian(p21,steps);
m22:=Smedian(p22,steps);
lowest:=min(min(m11,m12),min(m21,m22));
highest:=max(max(m11,m12),max(m21,m22));
result:=highest-lowest>100;
p11:=nil;
p12:=nil;
p21:=nil;
p22:=nil;
end;
function calc_weightF: double; {calculate weighting factor for different exposure duration and gain}
var
gain1,gain2 : double;
begin
if head.exposure<>0 then result:=head.exposure/head_ref.exposure else result:=1;{influence of each image depending on the exposure_time}
if head.egain<>head_ref.egain then {rare}
begin {check egain}
gain1:=strtofloat1(head_ref.egain);
gain2:=strtofloat1(head.egain);
if gain1<>0 then
result:=result*gain2/gain1; {-e/adu}
memo2_message('Warning different egain used!! '+copy(head.egain,1,5)+' ínstead of '+copy(head_ref.egain,1,5)+' [e-/ADU]. Will compensate accordingly.');
end
else
begin {check gain/iso}
if head.gain<>head_ref.gain then {rare}
memo2_message('Warning different gain used!! '+head.gain+' ínstead of '+head_ref.gain+'. Can not compensate unless egain [e-/ADU] is added manually to header.');
end;
end;
procedure stack_average(oversize:integer; var files_to_process : array of TfileToDo; out counter : integer);{stack average}
var
fitsX,fitsY,c,width_max, height_max,old_width, old_height,x_new,y_new,col,binning,oversizeV : integer;
background_correction, weightF,hfd_min : double;
init, solution,use_star_alignment,use_manual_align,use_ephemeris_alignment, use_astrometry_internal,vector_based : boolean;
tempval : single;
warning : string;
begin
with stackmenu1 do
begin
use_star_alignment:=stackmenu1.use_star_alignment1.checked;
use_manual_align:=stackmenu1.use_manual_alignment1.checked;
use_ephemeris_alignment:=stackmenu1.use_ephemeris_alignment1.checked;
use_astrometry_internal:=use_astrometry_internal1.checked;
hfd_min:=max(0.8 {two pixels},strtofloat2(stackmenu1.min_star_size_stacking1.caption){hfd});{to ignore hot pixels which are too small}
counter:=0;
sum_exp:=0;
jd_sum:=0;{sum of Julian midpoints}
jd_stop:=0;{end observations in Julian day}
init:=false;
background_correction:=0;
{simple average}
begin
for c:=0 to length(files_to_process)-1 do
if length(files_to_process[c].name)>0 then
begin
try { Do some lengthy operation }
ListView1.Selected :=nil; {remove any selection}
ListView1.ItemIndex := files_to_process[c].listviewindex;{show wich file is processed}
Listview1.Items[files_to_process[c].listviewindex].MakeVisible(False);{scroll to selected item}
filename2:=files_to_process[c].name;
Application.ProcessMessages;
{load image}
if ((esc_pressed) or (load_fits(filename2,true {light},true,init=false {update memo only for first ref img},0,head,img_loaded)=false)) then begin memo2_message('Error');{can't load} exit;end;{update memo for case esc is pressed}
if init=true then {first file done}
begin
if ((old_width<>head.width) or (old_height<>head.height)) then memo2_message('█ █ █ █ █ █ Warning different size image!');
if head.naxis3>length(img_average) {head.naxis3} then begin memo2_message('█ █ █ █ █ █ Abort!! Can'+#39+'t combine colour to mono files.'); exit;end;
end;
if init=false then
begin
binning:=report_binning(head.height);{select binning based on the height of the first light}
head_ref:=head;{backup solution}
initialise_var1;{set variables correct. Do this before apply dark}
initialise_var2;{set variables correct}
if ((bayerpat='') and (make_osc_color1.checked)) then
if stackmenu1.bayer_pattern1.Text='auto' then memo2_message('█ █ █ █ █ █ Warning, Bayer colour pattern not in the header! Check colours and if wrong set Bayer pattern manually in tab "stack alignment". █ █ █ █ █ █')
else
if test_bayer_matrix(img_loaded)=false then memo2_message('█ █ █ █ █ █ Warning, grayscale image converted to colour! Un-check option "convert OSC to colour". █ █ █ █ █ █');
end;
apply_dark_and_flat(img_loaded);{apply dark, flat if required, renew if different head.exposure or ccd temp}
{these global variables are passed-on in procedure to protect against overwriting}
memo2_message('Adding file: '+inttostr(c+1)+'-'+nr_selected1.caption+' "'+filename2+'" to average. Using '+inttostr(head.dark_count)+' darks, '+inttostr(head.flat_count)+' flats, '+inttostr(head.flatdark_count)+' flat-darks') ;
Application.ProcessMessages;
if esc_pressed then exit;
if make_osc_color1.checked then {do demosaic bayer}
begin
if head.naxis3>1 then memo2_message('█ █ █ █ █ █ Warning, light is already in colour ! Will skip demosaic. █ █ █ █ █ █')
else
demosaic_bayer(img_loaded); {convert OSC image to colour}
{head.naxis3 is now 3}
end;
if ((init=false ) and (use_astrometry_internal=false)) then {first image and not astrometry_internal}
begin
if ((use_manual_align) or (use_ephemeris_alignment)) then
begin
referenceX:=strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_X]); {reference offset}
referenceY:=strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_Y]); {reference offset}
end
else
begin
bin_and_find_stars(img_loaded, binning,1 {cropping},hfd_min,true{update hist},starlist1,warning);{bin, measure background, find stars}
find_quads(starlist1,0, quad_smallest,quad_star_distances1);{find quads for reference image}
pedestal_s:=cblack;{correct for difference in background, use cblack from first image as reference. Some images have very high background values up to 32000 with 6000 noise, so fixed pedestal_s of 1000 is not possible}
if pedestal_s<500 then pedestal_s:=500;{prevent image noise could go below zero}
background_correction:=pedestal_s-cblack;
head.datamax_org:=head.datamax_org+background_correction; if head.datamax_org>$FFFF then head.datamax_org:=$FFFF; {note head.datamax_org is already corrected in apply dark}
end;
end;
if init=false then {init}
begin
if oversize<0 then {shrink a lot, adapt in ratio}
begin
oversize:=max(oversize,-round((head.width-100)/2) );{minimum image width is 100}
oversizeV:=round(oversize*head.height/head.width);
height_max:=head.height+oversizeV*2;
end
else
begin
oversizeV:=oversize;
height_max:=head.height+oversize*2;
end;
width_max:=head.width+oversize*2;
setlength(img_average,head.naxis3,width_max,height_max);
setlength(img_temp,1,width_max,height_max);
for fitsY:=0 to height_max-1 do
for fitsX:=0 to width_max-1 do
begin
for col:=0 to head.naxis3-1 do
img_average[col,fitsX,fitsY]:=0; {clear img_average}
img_temp[0,fitsX,fitsY]:=0; {clear img_temp}
end;
old_width:=head.width;
old_height:=head.height;
if ((use_manual_align) or (use_ephemeris_alignment)) then
begin
referenceX:=strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_X]); {reference offset}
referenceY:=strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_Y]); {reference offset}
end;
end;{init, c=0}
solution:=true;
if use_astrometry_internal then sincos(head.dec0,SIN_dec0,COS_dec0) {do this in advance since it is for each pixel the same}
else
begin {align using star match}
if init=true then {second image}
begin
if ((use_manual_align) or (use_ephemeris_alignment)) then
begin {manual alignment}
solution_vectorX[2]:=referenceX-strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_X]); {calculate correction}
solution_vectorY[2]:=referenceY-strtofloat2(ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[L_Y]);
memo2_message('Solution x:=x+'+floattostr6(solution_vectorX[2])+', y:=y+'+floattostr6(solution_vectorY[2]));
end
else
begin{internal alignment}
bin_and_find_stars(img_loaded, binning,1 {cropping},hfd_min,true{update hist},starlist2,warning);{bin, measure background, find stars}
background_correction:=pedestal_s-cblack;
head.datamax_org:=head.datamax_org+background_correction; if head.datamax_org>$FFFF then head.datamax_org:=$FFFF; {note head.datamax_org is already corrected in apply dark}
find_quads(starlist2,0,quad_smallest,quad_star_distances2);{find star quads for new image}
if find_offset_and_rotation(3,strtofloat2(stackmenu1.quad_tolerance1.text)) then {find difference between ref image and new image}
memo2_message(inttostr(nr_references)+' of '+ inttostr(nr_references2)+' quads selected matching within '+stackmenu1.quad_tolerance1.text+' tolerance.'
+' Solution x:='+floattostr6(solution_vectorX[0])+'*x+ '+floattostr6(solution_vectorX[1])+'*y+ '+floattostr6(solution_vectorX[2])
+', y:='+floattostr6(solution_vectorY[0])+'*x+ '+floattostr6(solution_vectorY[1])+'*y+ '+floattostr6(solution_vectorY[2]) )
else
begin
memo2_message('Not enough quad matches <3 or inconsistent solution, skipping this image.');
files_to_process[c].name:=''; {remove file from list}
solution:=false;
ListView1.Items.item[files_to_process[c].listviewindex].SubitemImages[L_result]:=6;{mark 3th column with exclaimation}
ListView1.Items.item[files_to_process[c].listviewindex].subitems.Strings[2]:='no solution';{no stack result}
end;
end;{internal alignment}
end
else
reset_solution_vectors(1);{no influence on the first image}
end;
init:=true;{initialize for first image done}
if solution then
begin
inc(counter);
sum_exp:=sum_exp+head.exposure;
weightF:=calc_weightF;{calculate weighting factor for different exposure duration and gain}
date_to_jd(head.date_obs,head.exposure);{convert head.date_obs string and head.exposure time to global variables jd_start (julian day start head.exposure) and jd_mid (julian day middle of the head.exposure)}
if jd_start>jd_stop then jd_stop:=jd_start;{find latest start time}
jd_sum:=jd_sum+jd_mid;{sum julian days of images at midpoint head.exposure}
vector_based:=((use_star_alignment) or (use_manual_align) or (use_ephemeris_alignment));
if ((vector_based=false) and (a_order=0)) then {no SIP from astronomy.net}
begin
astrometric_to_vector;{convert astrometric solution to vector solution}
vector_based:=true;
end;
for fitsY:=1 to head.height do {skip outside "bad" pixels if mosaic mode}
for fitsX:=1 to head.width do
begin
calc_newx_newy(vector_based,fitsX,fitsY);{apply correction}
x_new_float:=x_new_float+oversize;y_new_float:=y_new_float+oversizeV;
x_new:=round(x_new_float);y_new:=round(y_new_float);
if ((x_new>=0) and (x_new<=width_max-1) and (y_new>=0) and (y_new<=height_max-1)) then
begin
for col:=0 to head.naxis3-1 do {all colors}
img_average[col,x_new,y_new]:=img_average[col,x_new,y_new]+ (img_loaded[col,fitsX-1,fitsY-1] +background_correction)*weightf;{image loaded is already corrected with dark and flat}{NOTE: fits count from 1, image from zero}
img_temp[0,x_new,y_new]:=img_temp[0,x_new,y_new]+weightF{typical 1};{count the number of image pixels added=samples.}
end;
end;
end;
progress_indicator(10+89*counter/(length(files_to_Process){ListView1.items.count}),' Stacking');{show progress}
finally
end;
end;
if counter<>0 then
begin
head_ref.naxis3:= head.naxis3; {store colour info in reference header}
head_ref.naxis:= head.naxis; {store colour info in reference header}
head:=head_ref;{restore solution variable of reference image for annotation and mount pointer. Works only if not resized}
head.height:=height_max;
head.width:=width_max;
setlength(img_loaded,head.naxis3,head.width,head.height);{new size}
For fitsY:=0 to head.height-1 do
for fitsX:=0 to head.width-1 do
begin {pixel loop}
tempval:=img_temp[0,fitsX,fitsY];
for col:=0 to head.naxis3-1 do
begin {colour loop}
if tempval<>0 then img_loaded[col,fitsX,fitsY]:=img_average[col,fitsX,fitsY]/tempval {scale to one image by diving by the number of pixels added}
else
begin { black spot filter or missing value filter due to image rotation}
if ((fitsX>0) and (img_temp[0,fitsX-1,fitsY]<>0)) then img_loaded[col,fitsX,fitsY]:=img_loaded[col,fitsX-1,fitsY]{take nearest pixel x-1 as replacement}
else
if ((fitsY>0) and (img_temp[0,fitsX,fitsY-1]<>0)) then img_loaded[col,fitsX,fitsY]:=img_loaded[col,fitsX,fitsY-1]{take nearest pixel y-1 as replacement}
else
img_loaded[col,fitsX,fitsY]:=0;{clear img_loaded since it is resized}
end; {black spot}
end;{colour loop}
end;{pixel loop}
end; {counter<>0}
// restore_solution(true);{restore solution variable of reference image for annotation and mount pointer}
end;{simple average}
end;{with stackmenu1}
{arrays will be nilled later. This is done for early exits}
end;
function minimum_distance_borders(fitsX,fitsY,w,h: integer): single;
begin
result:=min(fitsX,w-fitsX);
result:=min(fitsY,result);
result:=min(h-fitsY,result);
end;
procedure stack_mosaic(oversize:integer; var files_to_process : array of TfileToDo; max_dev_backgr: double; out counter : integer);{mosaic/tile mode}
var
fitsX,fitsY,c,width_max, height_max,x_new,y_new,col, cropW,cropH,iterations : integer;
value, dummy,median,median2,delta_median,correction,maxlevel,mean,noise : double;
tempval : single;
init, vector_based,merge_overlap,equalise_background : boolean;
background_correction,background_correction_center,background : array[0..2] of double;
counter_overlap : array[0..2] of integer;
begin
with stackmenu1 do
begin
{move often uses setting to booleans. Great speed improved if use in a loop and read many times}
merge_overlap:=merge_overlap1.checked;
Equalise_background:=Equalise_background1.checked;
counter:=0;
sum_exp:=0;
jd_sum:=0;{sum of Julian midpoints}
jd_stop:=0;{end observations in Julian day}
init:=false;
dummy:=0;
{mosaic mode}
begin
for c:=0 to length(files_to_process)-1 do
if length(files_to_process[c].name)>0 then
begin
try { Do some lengthy operation }
ListView1.Selected :=nil; {remove any selection}
ListView1.ItemIndex := files_to_process[c].listviewindex;{show wich file is processed}
Listview1.Items[files_to_process[c].listviewindex].MakeVisible(False);{scroll to selected item}
filename2:=files_to_process[c].name;
Application.ProcessMessages;
{load image}
if ((esc_pressed) or (load_fits(filename2,true {light},true,init=false {update memo only for first ref img},0,head,img_loaded)=false)) then begin memo2_message('Error');{can't load} exit;end;{update memo for case esc is pressed}
if init=true then
begin
// not for mosaic||| if init=true then if ((old_width<>head.width) or (old_height<>head.height)) then memo2_message('█ █ █ █ █ █ Warning different size image!');
if head.naxis3>length(img_average) {head.naxis3} then begin memo2_message('█ █ █ █ █ █ Abort!! Can'+#39+'t combine mono and colour files.'); exit;end;
end;
if init=false then
begin
head_ref:=head;{backup solution}
initialise_var1;{set variables correct}
initialise_var2;{set variables correct}
end;
memo2_message('Adding file: '+inttostr(c+1)+'-'+nr_selected1.caption+' "'+filename2+'" to mosaic.'); // Using '+inttostr(dark_count)+' dark(s), '+inttostr(flat_count)+' flat(s), '+inttostr(flatdark_count)+' flat-dark(s)') ;
Application.ProcessMessages;
if esc_pressed then exit;
if init=false then {init}
begin
oversize:=head.width*mosaic_width1.position div 2;{increase the oversize to have space for the tiles}
width_max:=head.width+oversize*2;
height_max:=head.height+oversize*2;
setlength(img_average,head.naxis3,width_max,height_max);
setlength(img_temp,1,width_max,height_max);{gray}
for fitsY:=0 to height_max-1 do
for fitsX:=0 to width_max-1 do
begin
for col:=0 to head.naxis3-1 do
begin
img_average[col,fitsX,fitsY]:=0; {clear img_average}
end;
img_temp[0,fitsX,fitsY]:=0; {clear img_temp}
end;
end;{init, c=0}
for col:=0 to head.naxis3-1 do {calculate background and noise if required}
begin
if equalise_background then
begin
background[col]:=mode(img_loaded,col,round(0.2*head.width),round(0.8*head.width),round(0.2*head.height),round(0.8*head.height),32000); {most common 80% center}
background_correction_center[col]:=1000 - background[col] ;
end
else
begin
background[col]:=0;
background_correction_center[col]:=0;
end;
end;
sincos(head.dec0,SIN_dec0,COS_dec0); {Alway astrometric. Do this in advance since it is for each pixel the same}
{solutions are already added in unit_stack}
begin
inc(counter);
sum_exp:=sum_exp+head.exposure;
date_to_jd(head.date_obs,head.exposure);{convert head.date_obs string and head.exposure time to global variables jd_start (julian day start head.exposure) and jd_mid (julian day middle of the head.exposure)}
if jd_start>jd_stop then jd_stop:=jd_start;{find latest start time}
jd_sum:=jd_sum+jd_mid;{sum julian days of images at midpoint head.exposure}
vector_based:=false;
if a_order=0 then {no SIP from astronomy.net}
begin
astrometric_to_vector;{convert astrometric solution to vector solution}
vector_based:=true;
end;
cropW:=trunc(stackmenu1.mosaic_crop1.Position*head.width/200);
cropH:=trunc(stackmenu1.mosaic_crop1.Position*head.height/200);
background_correction[0]:=0;
background_correction[1]:=0;
background_correction[2]:=0;
if init=true then {check image overlap intensisty differance}
begin
counter_overlap[0]:=0;
counter_overlap[1]:=0;
counter_overlap[2]:=0;
for fitsY:=1+1+cropH to head.height-1-cropH do {skip outside "bad" pixels if mosaic mode. Don't use the pixel at borders, so crop is minimum 1 pixel}
for fitsX:=1+1+cropW to head.width-1-cropW do
begin
calc_newx_newy(vector_based,fitsX,fitsY);{apply correction}
x_new:=round(x_new_float+oversize);y_new:=round(y_new_float+oversize);
if ((x_new>=0) and (x_new<=width_max-1) and (y_new>=0) and (y_new<=height_max-1)) then
begin
if img_loaded[0,fitsX-1,fitsY-1]>0.0001 then {not a black area around image}
begin
if img_average[0,x_new,y_new]<>0 then {filled pixel}
begin
for col:=0 to head.naxis3-1 do {all colors}
begin
correction:=round(img_average[col,x_new,y_new]-(img_loaded[col,fitsX-1,fitsY-1]+background_correction_center[col]) );
if abs(correction)<max_dev_backgr*1.5 then {acceptable offset based on the lowest and highest background measured earlier}
begin
background_correction[col]:=background_correction[col]+correction;
counter_overlap[col]:=counter_overlap[col]+1;
end;
end;
end;
end;
end;
end;
if counter_overlap[0]>0 then background_correction[0]:=background_correction[0]/counter_overlap[0];
if counter_overlap[1]>0 then background_correction[1]:=background_correction[1]/counter_overlap[1];
if counter_overlap[2]>0 then background_correction[2]:=background_correction[2]/counter_overlap[2];
end;
init:=true;{initialize for first image done}
for fitsY:=1+1+cropH to head.height-1-cropH do {skip outside "bad" pixels if mosaic mode. Don't use the pixel at borders, so crop is minimum 1 pixel}
for fitsX:=1+1+cropW to head.width-1-cropW do
begin
calc_newx_newy(vector_based,fitsX,fitsY);{apply correction}
x_new:=round(x_new_float+oversize);y_new:=round(y_new_float+oversize);
if ((x_new>=0) and (x_new<=width_max-1) and (y_new>=0) and (y_new<=height_max-1)) then
begin
if img_loaded[0,fitsX-1,fitsY-1]>0.0001 then {not a black area around image}
begin
dummy:=1+minimum_distance_borders(fitsX,fitsY,head.width,head.height);{minimum distance borders}
if img_temp[0,x_new,y_new]=0 then {blank pixel}
begin
for col:=0 to head.naxis3-1 do {all colors}
img_average[col,x_new,y_new]:=img_loaded[col,fitsX-1,fitsY-1]+background_correction_center[col] +background_correction[col];{image loaded is already corrected with dark and flat}{NOTE: fits count from 1, image from zero}
img_temp[0,x_new,y_new]:=dummy;
end
else
begin {already pixel filled, try to make an average}
for col:=0 to head.naxis3-1 do {all colors}
begin
median:=background_correction_center[col] +background_correction[col]+median_background(img_loaded,col,15,15,fitsX-1,fitsY-1);{find median value in sizeXsize matrix of img_loaded}
if merge_overlap=false then {method 2}
begin
median2:=median_background(img_average,col,15,15,x_new,y_new);{find median value of the destignation img_average}
delta_median:=median-median2;
img_average[col,x_new,y_new]:= img_average[col,x_new,y_new]+ delta_median*(1-img_temp[0,x_new,y_new]{distance border}/(dummy+img_temp[0,x_new,y_new]));{adapt overlap}
end
else
begin {method 1}
value:=img_loaded[col,fitsX-1,fitsY-1]+background_correction_center[col];
local_sd(fitsX-1-15 ,fitsY-1-15, fitsX-1+15,fitsY-1+15,col,img_loaded, {var} noise,mean, iterations);{local noise recheck every 10 th pixel}
maxlevel:=median+noise*5;