-
Notifications
You must be signed in to change notification settings - Fork 9
/
Copy pathmirone.m
5854 lines (5428 loc) · 280 KB
/
mirone.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
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
function varargout = mirone(varargin)
% MIRONE, by itself, creates a window bar from which you load a lot of grid/images formats
% MIRONE(FNAME) opens the file FNAME and displays it
% H = MIRONE(...) returns the handle to a new mirone window
%
% mirone('CALLBACK',handles,...) calls the local function named CALLBACK with the given input arguments.
% Copyright (c) 2004-2022 by J. Luis
%
% This program is part of Mirone and is free software; you can redistribute
% it and/or modify it under the terms of the GNU Lesser General Public
% License as published by the Free Software Foundation; either
% version 2.1 of the License, or any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
% Lesser General Public License for more details.
%
% Contact info: w3.ualg.pt/~jluis/mirone
% --------------------------------------------------------------------
% $Id: mirone.m 11437 2019-07-18 15:20:44Z j $
if (nargin > 1 && ischar(varargin{1}))
if (~isempty(strfind(varargin{1},':')) || ~isempty(strfind(varargin{1},filesep)) ) % A file name
if (ischar(varargin{2}))
% OK, this may mean a name with those horrendous stupid blanks. And we don't know how many
% So scan the first n elemennts on nargin that are chars, cat the and try to see if they much to
% an existing file. If it does, assume that's fname in, otherwise we assume that the remaing
% args are of the form -C<...> -X<...>, which are allowed for calling the external_drive() mechanism.
n = 2; % At least this we know
for (k = 3:nargin)
if (~ischar(varargin{k})), break, end % Increments conter on number of first args being char
n = n + 1;
end
for (k = 2:n)
fmt = repmat('%s ', 1, k); fmt(end) = []; % We don't want the last ' ' in name
t = sprintf(fmt, varargin{1:k});
if (exist(t, 'file'))
varargin{1} = t; varargin(2:k) = [];
break
end
end
end
h = mirone_OpeningFcn(varargin{:});
if (nargout), varargout{1} = h; end
else
if (varargin{1}(1) ~= '-') % Make sure the argument is not any of -C..., -X... magic args
gui_Callback = str2func(varargin{1});
[varargout{1:nargout}] = feval(gui_Callback,varargin{2:end});
else
[varargout{1:nargout}] = mirone_OpeningFcn(varargin{:});
end
end
else
h = mirone_OpeningFcn(varargin{:});
if (nargout), varargout{1} = h; end
end
% --------------------------------------------------------------------------------------------------
function hObject = mirone_OpeningFcn(varargin)
% PRAGMA SECTION (It's far far from clear when files must be declared here)
%#function uigetfolder_standalone coordinate_system gmtmbgrid_m gmtedit
%#function grdtrack_m
%#function grdutils scaleto8 bpass3d inv3d rtp3d syn3d igrf_m
%#function range_change swan tsun2 mansinha_m deform_mansinha deform_okada dim_funs
%----- These are for image
%#function grayto8 grayto16 grayxform imfilter_mex imhistc imlincombc parityscan intlutc ordf
%#function imreconstructmex applylutc bwboundariesmex bwlabel1 bwlabel2 ditherc bwlabelnmex morphmex
%----- These are in utils
%#function tabpanelfcn degree2dms dms2degree dec2deg dec_year ivan_the_terrible ddewhite string_token
%#function test_dms text_read double2ascii save_seismicity jd2date guess_file shading_mat
%#function getline_j new_frame3D histos_seis uirestore_j uisuspend_j
%----- These is for write_gmt_script
%#function draw_scale time_stamp pscoast_options_Mir paint_option w_option
%----- Those are ..., the hell with explanations for what I don't realy understand. They are needed, that's all.
%#function country_select read_isf choosebox magbarcode listbox_message add_poles animate_seismicity
%#function get_polygon rot_euler datums telha_m find_clusters fft_stuff select_cols uistack_j smoothing_param
%#function patch_meca ui_edit_patch_special bands_list multibandread_j imscroll_j geog2projected_pts
%#function mltable_j iptcheckinput resampsep wgifc telhometro vitrinite edit_line move_obj make_arrow
%#function edit_track_mb save_track_mb houghmex qhullmx writegif mpgwrite cq helpdlg
%#function move2side aguentabar gdal_project gdalwarp_mex poly2mask_fig url2image calc_bonin_euler_pole spline_interp
%#function mat2clip buffer_j PolygonClip akimaspline shake_mex ground_motion wms_tool microlev
%#function write_esri_hdr distmin mag_synthetic image_histo write_gmt_symb mkpj decompress mosaicer
%#function lasreader_mex laszreader_mex escorrega show_manguito travel thresholdit intersections nswing runCB_tintol
%#function usgs_recent_seismicity sat_orbits uisetdate doy sshist win_open_mex show_mb mb_cleaning_params
%#function c_cpt2cmap c_grdfilter c_grdinfo c_grdlandmask c_grdproject c_grdread c_grdsample
%#function c_grdtrend c_mapproject c_nearneighbor c_shoredump c_surface popenr diffCenterVar hellinger bingham
%#function gmtlist_m grdgradient_m trend1d_m external_drive chimoce interp_chimoce earth_tides earthtide DateStr2Num
%#function ddist eucdist2 masker nth_element filter1d cayula_cornillon
global home_dir; fsep = filesep;
toCompile = false; % To compile set this one to TRUE
if (toCompile)
home_dir = cd;
else
if (isempty(home_dir)) % First time call. Find out where we are
home_dir = fileparts(mfilename('fullpath')); % Get the Mirone home dir and set path
addpath(home_dir, [home_dir fsep 'src_figs'],[home_dir fsep 'utils']);
if (exist('OCTAVE_VERSION','builtin') ~= 0) % This is a repetition of the test later in mirone_uis
addpath([home_dir fsep 'lib_mex' fsep 'octave' fsep octave_config_info.canonical_host_type]);
else
addpath([home_dir fsep 'lib_mex']);
end
end
end
[hObject,handles,home_dir] = mirone_uis(home_dir);
pos = get(handles.figure1,'pos');
if (pos(4) == 32) % On Win 8 and retina, the bloody OS insists in resetting it to 32
set(handles.figure1,'Pos',[pos(1), pos(2)+31, pos(3), 1]);
end
handles.home_dir = home_dir;
handles.DefLineThick = 0.5; % Default line thickness (overwriten by mirone_pref)
handles.DefLineColor = 'k'; % Default line color (overwriten by mirone_pref)
handles.DefineMeasureUnit = 'k';% Default measure units to kilometrs (overwriten by mirone_pref)
handles.grdname = []; % Contains the name of the current (if it's the case) gmt grid
handles.nTrack = 0; % Counter of the number of MB tracks imported
handles.origFig = []; % To store the original image copy
handles.fileName = []; % To store any input grid/image file name
handles.image_type = 0; % 1->grd; 2-> trivial (jpg,png,bmp,etc...); 3->GeoTIFF; 4->DEMs; 20-> white bg
handles.transparency = []; % Will hold the word 'alpha' when reading a RGBA image
handles.computed_grid = 0; % However, matrices with a gmt header will have this == 1, so that they can be saved
handles.no_file = 1; % 0 means a grid is loaded and 1 that it is not (to test when icons can be pushed)
handles.geog = 1; % By default grids are assumed to be in geographical coordinates
handles.swathRatio = 3; % Default swath width / water depth ratio for multibeam planing
handles.grdMaxSize = 524288000; % I use this for limiting the grid size that is stored in RAM (50 Mb)
handles.firstMBtrack = 1; % Used for knowing whether to display or not the MB planing info message in "start planing"
handles.EarthRad = 6371.005076; % Authalic radius
handles.maregraphs_count = 0; % Counter of maregraphs (tsunami modeling)
handles.Illumin_type = 0; % 0 -> no illumination; 1 -> grdgradient; 2 -> grdgradient Lambertian; 4 -> Lambertian;
handles.zoom_state = 0; % Flag to signal if zoom state is to be re-activated (0 means no re-activation)
handles.bg_color = [1 1 1]; % (default) Backgoround color used when grid values == NaN
handles.which_cont = []; % To store the contour levels (used to not draw repeated contours)
handles.have_nans = 0; % Used to know if the grids have NaNs
handles.is_draped = false; % Used to know if the image comes from draping
handles.was_int16 = false; % Keep track of short int imported grids
handles.was_pixreg = false; % Keep track of the grid registration type. Currently, only used in plot_composer
handles.Nodata_int16 = []; % To store Nodata of short int grids
handles.ForceInsitu = false; % Use "insitu" grid importing (transposition)
handles.DefineEllipsoide = [6378137, 0, 1/298.2572235630]; % Defaults to WGS-84
handles.path_tmp = [home_dir fsep 'tmp' fsep];
handles.path_data = [home_dir fsep 'data' fsep];
handles.last_directories = {handles.path_tmp; home_dir}; % Let it have something existent
handles.hImg = []; % To hold the image handle
handles.firstIllum = 1; % First illumination will use the displayed image which may have been IP
handles.flederBurn = 1; % When build a fleder obj burn eventual coastlines
handles.flederPlanar = 1; % Default to planar flder objs (but mirone_pref knows better)
handles.whichFleder = 1; % whichFleder = 1 for the free iview4d or 0 for the true thing (fledermaus)
handles.oldSize = [get(hObject,'Pos'); get(hObject,'Pos')]; % Duplicate so that we store ORIGINAL size
if (handles.oldSize(1,4) == 0), handles.oldSize(1,4) = 1; end
handles.is_projected = false; % To keep track if coords are projected or not
handles.defCoordsIn = 0; % To use when Load files and have to decide if we need to project
% 0 -> don't know; -1 -> Coords are already projected; 1 -> geog coords needing project
try handles.Projections; % Use a try/catch since isfield is brain-dead long
catch, handles.Projections = nan; % Set it to something to prevent "unknown field" error
end
handles.scale2meanLat = 1; % Scale geog ref images so that at middle image dx = dy in cartesian units (kind of proj)
handles.FOpenList = cell(numel(handles.RecentF),1);
handles.withSliders = true; % Set Zoom sliders
handles.validGrid = false; %
handles.nLayers = 1; % If > 1 after reading a netCDF file call aquamoto
handles.deflation_level = 0; % If > 0 will create compressed netCDF-4 files
handles.is_defRegion = false; % A def region is a particular case to create GMT custom symbols.
handles.screenSize = []; % Fill it for the case when one wants to fake a different get(0, 'ScreenSize')
handles.img_with_minmax = []; % Will have [lower higher] when scaleto8 is to operate on an interval
handles.cmapIsLog = 0; % So that we know how to update logaritmized cmaps.
try % A file named mirone_pref.mat contains the preferences, read them from it
prf = load([handles.path_data 'mirone_pref.mat']);
handles.geog = prf.geog;
handles.grdMaxSize = prf.grdMaxSize; % 2^20 = 1 Mb
handles.swathRatio = prf.swathRatio;
handles.last_directories = prf.directory_list;
handles.DefLineThick = sscanf(strtok(prf.DefLineThick{1}),'%f');
% Decode the line color string into the corresponding char (e.g. k,w, etc...)
if (strcmp(prf.DefLineColor{1},'Black')), handles.DefLineColor = 'k';
else, handles.DefLineColor = lower(prf.DefLineColor{1}(1));
end
% Decode the Measure unities into a char code (e.g. n, k, m, u from {'nautical miles' 'kilometers' 'meters' 'user'})
handles.DefineMeasureUnit = prf.DefineMeasureUnit{1}(1);
handles.DefineEllipsoide = prf.DefineEllipsoide_params; % Set the default ellipsoide parameters (a,b,f)
handles.flederBurn = prf.flederBurn;
handles.flederPlanar = prf.flederPlanar;
handles.scale2meanLat = prf.scale2meanLat;
handles.FOpenList = prf.FOpenList;
handles.whichFleder = prf.whichFleder;
if (~prf.moveDoubleClick) % this info is used by UI_EDIT_POLYGON()
setappdata(handles.axes1,'MovPolyg','extend') % Move lines with a Shift-click drag-n-drop
end
handles.bg_color = prf.nanColor;
handles.deflation_level = prf.deflation_level;
if (isfield(prf, 'TDRver')), handles.TDRver = prf.TDRver;
else, handles.TDRver = '2.0';
end
catch
% Tell mirone_pref to write up the defaults.
mirone_pref(handles,'nikles')
handles = guidata(handles.figure1); % And need also the updated handles
end
j = false(1,numel(handles.last_directories)); % vector for eventual cleaning non-existing dirs
for (i = 1:numel(handles.last_directories)) % Check that all dirs in last_directories exist
j(i) = (exist(handles.last_directories{i},'dir') ~= 7);
end
handles.last_directories(j) = []; % clean non-existing directories
if (isempty(handles.last_directories)) % Don't ever let it be empty
handles.last_directories = {handles.path_tmp; home_dir}; % Let it have something existent
end
if (any(j)) % If any of the old dirs evaporated, update that info in mirone_prefs
directory_list = handles.last_directories; isempty(directory_list); % To shut up MLint
if (handles.version7), save([handles.path_data 'mirone_pref.mat'],'directory_list', '-append', '-v6')
else, save([handles.path_data 'mirone_pref.mat'],'directory_list', '-append')
end
end
handles.work_dir = handles.last_directories{1};
if (handles.work_dir ~= fsep), handles.work_dir = [handles.work_dir fsep]; end
handles.last_dir = handles.last_directories{1}; % Initialize last_dir to work_dir
setappdata(hObject,'swathRatio',handles.swathRatio); % I need this for getline_mb
% Change the MeasureDistance label to the selected (in prefs) units
set(handles.ToolsMeasureDist,'Label',['Distance in ' handles.DefineMeasureUnit])
% -------------------------- Find in which mode Mirone was called ----------------------------------
drv = []; grd_data_in = 0; grd_data_interfero = 0; pal = []; win_name = 'Mirone';
feval_me = ''; x_comm = cell(1,numel(varargin));
% Pre-process vargin to detect, and extract, "self-driving" commands in the form -C<...>
if (~isempty(varargin))
c = false(numel(varargin),1);
for (k = 1:numel(varargin))
if (ischar(varargin{k}) && varargin{k}(1) == '-')
if (varargin{k}(2) == 'C')
feval_me = varargin{k}(3:end); % A -C<...> command to be executed by feval
else
x_comm{k} = varargin{k}(3:end); % A -X<...> command to be send to child
end
c(k) = true;
end
end
if (any(c)), varargin(c) = []; end
end
if (~isempty(varargin))
n_argin = numel(varargin);
if (n_argin == 1 && ischar(varargin{1})) % Called with a file name as argument
[pato, fname, EXT] = fileparts(varargin{1}); % Test to check online command input
if (isempty(pato))
varargin{1} = [cd fsep fname EXT];
handles.last_dir = cd;
else
handles.last_dir = pato;
end
[drv, algures] = aux_funs('findFileType',varargin{1});
if (ischar(algures)), varargin{1} = algures; end % File exists but not in Mirone's root dir
if (ischar(algures) || algures), handles.fileName = varargin{1}; end % Can be added by recentFiles
elseif (isa(varargin{1},'uint8') || isa(varargin{1},'int8') || islogical(varargin{1}))
% Called with an image as argument and optionaly an struct header (& geog, name, cmap optional fields)
if (isa(varargin{1},'int8')) % We cannot represent a int8 image. Do something
varargin{1} = uint8(cvlib_mex('addS', int16(varargin{1}), 128)); % [-128 127] -> [0 255]
end
% Now deal with the case of a eventual multiband ( > than 3 planes) array
if (size(varargin{1},3) > 3), aux_funs('toBandsList', handles.figure1, varargin{1}, 'multiband array'), end
isReferenced = false;
if (n_argin >= 2 && isa(varargin{2},'struct')) % An image with coordinates
tmp = varargin{2};
handles.head = tmp.head; X = tmp.X; Y = tmp.Y;
handles.image_type = 3; axis_t = 'xy';
if (isfield(tmp,'geog')), handles.geog = tmp.geog; end % Prevails over the guess in show_image (does it?)
if (isfield(tmp,'cmap')), set(handles.figure1,'Colormap',tmp.cmap); end
if (isfield(tmp,'name')), win_name = tmp.name; end
if (isfield(tmp,'srsWKT'))
aux_funs('appP', handles, tmp.srsWKT) % If we have a WKT proj, store it
isReferenced = true;
if (~handles.geog), handles.is_projected = true; end % WEAK LOGIC. SHOULD PARSE WKT TO MAKE SURE
else
projWKT ='';
if (n_argin == 3 && ishandle(varargin{3}))
handTmp = guidata(varargin{3});
projWKT = getappdata(handTmp.figure1,'ProjWKT');
if (isempty(projWKT)), projWKT = getappdata(handTmp.figure1,'Proj4'); end
handles.geog = handTmp.geog + 10;
if (isappdata(handTmp.figure1,'GCPregImage'))
setappdata(hObject,'GCPregImage',getappdata(handTmp.figure1,'GCPregImage'))
setappdata(hObject,'fnameGCP', win_name)
end
end
if (~isempty(projWKT)), isReferenced = true; tmp.srsWKT = projWKT; end
end
else
X = []; Y = []; win_name = 'Cropped_image';
handles.image_type = 2; handles.geog = 0; axis_t = 'off';
handles.head = [1 size(varargin{1}, 2) 1 size(varargin{1}, 1) 0 255 0 1 1]; % Fake a grid reg GMT header
if (ndims(varargin{1}) == 2), set(handles.figure1,'Colormap',gray(256)); end
pal = getappdata(0,'CropedColormap'); % See if we have a colormap to use here
if (~isempty(pal)), set(handles.figure1,'Colormap',pal), rmappdata(0,'CropedColormap'), end
setappdata(hObject,'Croped','yes'); % ???
end
if (size(varargin{1},3) == 4), handles.transparency = 'alpha'; end % Signal show_image that 4th band is transparency
handles = show_image(handles,win_name,X,Y,varargin{1},0,axis_t,handles.head(7),1);
if (~isReferenced), grid_info(handles,[],'iminfo',varargin{1}); % Create a info string
else, grid_info(handles,tmp.srsWKT,'referenced',varargin{1});
end
handles = aux_funs('isProj',handles); % Check/set about coordinates type
elseif (n_argin == 1 && isa(varargin{1},'struct') && isfield(varargin{1},'proj4'))
% A GMT grid/image structure. (for images we still do not use eventual alpha channel)
handles.head = [varargin{1}.range varargin{1}.registration varargin{1}.inc];
win_name = 'Nikles';
if (~isempty(varargin{1}.title)), win_name = varargin{1}.title; end
if (~isfield(varargin{1}, 'image'))
Z = varargin{1}.z; grd_data_in = true;
if (~isa(Z,'single')), Z = single(Z); end
handles.have_nans = grdutils(Z,'-N');
X = linspace(varargin{1}.range(1), varargin{1}.range(2), size(varargin{1}.z, 2));
Y = linspace(varargin{1}.range(3), varargin{1}.range(4), size(varargin{1}.z, 1));
else
handles.image_type = 2; isReferenced = false;
X = [varargin{1}.x(1) varargin{1}.x(end)]; Y = [varargin{1}.y(1) varargin{1}.y(end)];
handles.geog = aux_funs('guessGeog', [Y Y]);
ProjectionRefWKT = varargin{1}.wkt;
if (isempty(ProjectionRefWKT) && ~isempty(varargin{1}.proj4))
ProjectionRefWKT = ogrproj(varargin{1}.proj4);
end
if (~isempty(ProjectionRefWKT))
aux_funs('appP', handles, varargin{1}.wkt) % If we have a WKT proj, store it
isReferenced = true;
if (~handles.geog), handles.is_projected = true; end % WEAK LOGIC. SHOULD PARSE WKT TO MAKE SURE
end
ax_t = 'off';
if (handles.geog || isReferenced)
handles.image_type = 3; ax_t = 'xy';
end
if (ndims(varargin{1}.image) == 2), set(handles.figure1,'Colormap',gray(256)); end
alpha = 1;
if (~isempty(varargin{1}.alpha)), alpha = varargin{1}.alpha; end
handles = show_image(handles,win_name,X,Y,varargin{1}.image,0,ax_t,varargin{1}.registration,1, alpha);
if (~isReferenced), grid_info(handles,[],'iminfo',varargin{1}.image); % Create a info string
else, grid_info(handles, ProjectionRefWKT, 'referenced', varargin{1}.image);
end
handles = aux_funs('isProj',handles); % Check/set about coordinates type
end
elseif (n_argin == 2 && isequal([size(varargin{2},1) size(varargin{2},2)],[1 9]))
% A matrix with the classic nine elements header vector
handles.head = varargin{2};
Z = varargin{1}; grd_data_in = true;
if (~isa(Z,'single')), Z = single(Z); end
handles.have_nans = grdutils(Z,'-N');
X = linspace(handles.head(1), handles.head(2), size(Z,2));
Y = linspace(handles.head(3), handles.head(4), size(Z,1));
elseif (n_argin == 3 && numel(varargin{1}) == size(varargin{3},2) && numel(varargin{2}) == size(varargin{3},1))
% A (X, Y, Z) where X,Y are the coord vectors and Z is the grid matrix
Z = varargin{3}; grd_data_in = true;
X = varargin{1}(:)'; Y = varargin{2}(:)';
if (X(2) < X(1)), X = X(end:-1:1); end
if (Y(2) < Y(1)), Y = Y(end:-1:1); end
if (~isa(Z,'single')), Z = single(Z); end
handles.have_nans = grdutils(Z,'-N');
win_name = 'Nikles';
zz = grdutils(Z,'-L');
handles.head = [X(1) X(end) Y(1) Y(end) zz(1) zz(2) 0 X(2)-X(1) Y(2)-Y(1)];
elseif (n_argin < 4 && ~(isa(varargin{1},'uint8') || isa(varargin{1},'int8')))
% A matrix with an optional header. Treat it as a gmt grid. No error testing on the grid head descriptor
% A third argument is optional but if exists must contain a Mirone's fig or axes handle.
% The prefered way of calling this option should now be:
% mirone(Z, struct('head',[...], 'name','...'), [parent_mir_fig_handle])
Z = varargin{1}; grd_data_in = true;
if (~isa(Z,'single') && ~isa(Z,'uint16') && ~isa(Z,'int16')), Z = single(Z); end
handles.have_nans = grdutils(Z,'-N');
if (numel(varargin) >= 2 && isa(varargin{2},'struct')) % A grid with a header
tmp = varargin{2};
handles.head = tmp.head;
if (~isfield(tmp,'X'))
X = linspace(handles.head(1), handles.head(2), size(Z,2));
Y = linspace(handles.head(3), handles.head(4), size(Z,1));
else
X = tmp.X; Y = tmp.Y;
end
if (isfield(tmp,'name')), win_name = tmp.name; end % All calls should transmit a name, but ...
projWKT ='';
if (isfield(tmp,'srsWKT')), projWKT = tmp.srsWKT; end % If exists, take precedence
if (n_argin == 3 && ishandle(varargin{3}))
handTmp = guidata(varargin{3});
pal = get(handTmp.figure1, 'colormap');
if (isempty(projWKT)), projWKT = getappdata(handTmp.figure1,'ProjWKT'); end
if (isempty(projWKT)), projWKT = getappdata(handTmp.figure1,'Proj4'); end
handles.geog = handTmp.geog + 10;
if (isappdata(handTmp.figure1,'GCPregImage'))
setappdata(hObject,'GCPregImage',getappdata(handTmp.figure1,'GCPregImage'))
setappdata(hObject,'fnameGCP', win_name)
end
end
if (isfield(tmp,'geog')), handles.geog = tmp.geog + 10; end % If exists, take precedence
if (~isempty(projWKT))
grid_info(handles,projWKT,'referenced',varargin{1}); % Create a info string and save projWKT
handles.is_projected = true; % WEAK LOGIC. SHOULD PARSE WKT TO MAKE SURE
elseif (isfield(tmp,'ProjGMT')) % From geog_calculator. Has opt_J.
projection_menu(handles, tmp.ProjGMT)
handles = guidata(hObject); % Get the updated version changed in the above call
end
if (isfield(tmp,'cmap')), pal = tmp.cmap; end % Takes precedence over the above.
if (isfield(tmp,'screenSize'))
handles.screenSize = tmp.screenSize;
end
else
zz = grdutils(Z,'-L');
handles.head = [1 size(Z,2) 1 size(Z,1) zz(1) zz(2) 0 1 1];
X = 1:size(Z,2); Y = 1:size(Z,1);
end
elseif (n_argin == 4 && isnumeric(varargin{1}) && isa(varargin{2},'struct') && ...
(strcmp(varargin{3},'Deformation') || strcmp(varargin{3},'Interfero')) )
% A matrix. Treat it as if it'is a gmt grid. No error testing on the grid head descriptor
% Note: this is a special case of the situation above that will be used to identify this figure
% as an Okada deformtion data (via its Name). This info is searched by the tsunami modeling option
if (~isa(varargin{1},'single')), varargin{1} = single(varargin{1}); end
handles.head = varargin{2}.head; X = varargin{2}.X; Y = varargin{2}.Y;
Z = varargin{1};
handles.have_nans = grdutils(Z,'-N');
if (varargin{3}(1) == 'D')
grd_data_in = true; win_name = 'Okada deformation';
setappdata(hObject,'hFigParent',varargin{4});
else % A matrix input containing an interfeogram with cdo == varargin{4}
grd_data_interfero = true; win_name = 'Interferogram';
cdo = varargin{4};
end
end
end
% The following IF cases deal only with cases where a grid was given in argument
if (grd_data_in || grd_data_interfero)
handles.image_type = 1; handles.computed_grid = 1; % Signal that this is a computed grid
if (isempty(pal)), pal = jet(256); end
if (grd_data_interfero) % Interferogram grid
pal = load([handles.path_data 'gmt_other_palettes.mat'],'circular'); pal = pal.circular;
zz = uint8(abs(rem(double(Z),cdo)/cdo)*255);
else
zz = scaleto8(Z);
end
set(handles.figure1,'Colormap',pal)
aux_funs('StoreZ',handles,X,Y,Z) % If grid size is not to big we'll store it
aux_funs('colormap_bg',handles,Z,pal);
handles = show_image(handles,win_name,X,Y,zz,1,'xy',handles.head(7));
end
% ---------+================+---------- FIGURE VISIBLE HERE ----------+==============+-----------
set(handles.figure1,'Vis', 'on')
% ---------+================+---------- FIGURE VISIBLE HERE ----------+==============+-----------
handles.IAmAMac = strncmp(computer,'MAC',3);
setappdata(0,'IAmAMac',handles.IAmAMac)
if (handles.IAmAMac)
if (isempty(getappdata(0,'have_DYLD'))) % Deal with MacOSX blindness
%set_gmt(['DYLD_LIBRARY_PATH=' cd ':'],'DYLD_LIBRARY_PATH') % Prepend current dir to DYLD_LIBRARY_PATH
DYLD = getenv('DYLD_LIBRARY_PATH'); % Available only on > R14?? versions
setenv('DYLD_LIBRARY_PATH',[DYLD ':' cd])
setappdata(0,'have_DYLD',true) % Signal to do this only once
end
% F. TMW just cannot make things work decently on Macs. Labels are bigger than reserved space
% And now they screw it even more. The following simply f... on 2012 (and maybe earlier(
v = version('-release');
if (str2double(v(1:4)) < 2012)
set(handles.File, 'Label', ['<HTML><FONT size="3">' get(handles.File, 'Label') '</Font></html>'])
set(handles.Image, 'Label', ['<HTML><FONT size="3">' get(handles.Image, 'Label') '</Font></html>'])
set(handles.Tools, 'Label', ['<HTML><FONT size="3">' get(handles.Tools, 'Label') '</Font></html>'])
set(handles.Draw, 'Label', ['<HTML><FONT size="3">' get(handles.Draw, 'Label') '</Font></html>'])
set(handles.Geography, 'Label', ['<HTML><FONT size="3">' get(handles.Geography, 'Label') '</Font></html>'])
set(handles.Plates,'Label', ['<HTML><FONT size="3">' get(handles.Plates, 'Label') '</Font></html>'])
set(handles.MagGrav, 'Label', ['<HTML><FONT size="3">' get(handles.MagGrav, 'Label') '</Font></html>'])
set(handles.Seismology, 'Label', ['<HTML><FONT size="3">' get(handles.Seismology, 'Label') '</Font></html>'])
set(handles.Tsunamis, 'Label', ['<HTML><FONT size="3">' get(handles.Tsunamis, 'Label') '</Font></html>'])
set(handles.GMT, 'Label', ['<HTML><FONT size="3">' get(handles.GMT, 'Label') '</Font></html>'])
set(handles.MetOce, 'Label', ['<HTML><FONT size="3">' get(handles.MetOce, 'Label') '</Font></html>'])
set(handles.GridTools, 'Label', ['<HTML><FONT size="3">' get(handles.GridTools, 'Label') '</Font></html>'])
set(handles.Projections, 'Label', ['<HTML><FONT size="3">' get(handles.Projections, 'Label') '</Font></html>'])
set(handles.Help, 'Label', ['<HTML><FONT size="3">' get(handles.Help, 'Label') '</Font></html>'])
end
set(hObject,'DockControls','off')
end
%Find out which gmt version is beeing used.
info = getappdata(0,'gmt_version'); % See if the info is already here.
if (isempty(info))
pato = [home_dir fsep 'gmt_share/coast/'];
info.full = (exist([pato 'binned_GSHHS_f.nc'], 'file') == 2);
info.high = (exist([pato 'binned_GSHHS_h.nc'], 'file') == 2);
info.intermediate = (exist([pato 'binned_GSHHS_i.nc'], 'file') == 2);
info.low = (exist([pato 'binned_GSHHS_l.nc'], 'file') == 2);
setappdata(0,'gmt_version',info); % Save it to the next time a new mirone window is opened
if (info.intermediate && info.low && (exist([home_dir fsep 'gmt_share/cpt/'], 'dir') == 7))
set_gmt(['GMT5_SHAREDIR=' home_dir fsep 'gmt_share']); % GMT5_SHAREDIR because it's searched before GMT_SHAREDIR
end
end
if (~info.full && ~info.high && ~info.intermediate && ~info.low)
info = set_gmt(['GMT_USERDIR=' home_dir fsep 'gmt_userdir']); % Even this one will have to go
setappdata(0,'gmt_version',info); % Save it to the next time a new mirone window is opened
end
if (~info.full && ~strcmp(info.full, 'y'))
set([handles.CoastLineFull handles.PBFull handles.RiversFull], 'Enable','off')
end
if (~info.high && ~strcmp(info.high, 'y'))
set([handles.CoastLineHigh handles.PBHigh handles.RiversHigh], 'Enable','off')
end
if (~info.intermediate && ~strcmp(info.intermediate, 'y'))
set([handles.CoastLineInterm handles.PBInterm handles.RiversInterm], 'Enable','off')
end
if (~info.low && ~strcmp(info.low, 'y'))
set([handles.CoastLineLow handles.PBLow handles.RiversLow], 'Enable','off')
end
% Deal with the new troubles introduced by using GMT5.2 that needs to know where to find its own share dir
t = set_gmt('GMT5_SHAREDIR', 'whatever'); % Inquire if GMT5_SHAREDIR exists
if (isempty(t)), t = set_gmt('GMT_SHAREDIR', 'whatever'); end % Then GMT_SHAREDIR
if (isempty(t) || (~(exist(t, 'dir') == 7))) % Test also that the dir exists
% If not, set a fake one with minimalist files so that GMT does not complain/errors
% We have to use GMT5_SHAREDIR and NOT GMT_SHAREDIR because it's the first one checked in gmt_init.c/GMT_set_env()
set_gmt(['GMT5_SHAREDIR=' home_dir fsep 'gmt_share']);
end
% See if we have a TMP dir set by a ENV var that will take precedence over the default one
t = set_gmt('MIRONE_TMP', 'whatever'); % Inquire if MIRONE_TMP exists
if (~isempty(t)) % Now check that the dir exists
if (exist(t, 'dir') == 7)
handles.path_tmp = [t fsep];
end
end
set_gmt(['PROJ_LIB=' home_dir fsep 'gmt_share' fsep 'proj']); % For projections with GDAL
set_gmt(['GDAL_DATA=' home_dir fsep 'gmt_share' fsep 'GDAL_DATA']);
%set_gmt(['PROJ_LIB=' home_dir fsep 'data' fsep 'proj_lib']); % For projections with GDAL
%set_gmt(['GDAL_DATA=' home_dir fsep 'data' fsep 'gdal_data']);
%set_gmt(['GEOTIFF_CSV=' home_dir fsep 'data' fsep 'gdal_data']);
set_gmt(['GEOTIFF_CSV=' home_dir fsep 'gmt_share' fsep 'GDAL_DATA']);
set_gmt('GDAL_HTTP_UNSAFESSL=YES'); % To tell libcurl to not verify the peer
%set_gmt('PROJ_DEBUG=5');
%set_gmt('CPL_DEBUG=ON');
guidata(hObject, handles);
t.home_dir = home_dir; t.work_dir = handles.work_dir; t.last_dir = handles.last_dir; t.version7 = handles.version7;
t.path_data = handles.path_data; t.path_tmp = handles.path_tmp; t.IamCompiled = handles.IamCompiled;
setappdata(0,'MIRONE_DIRS',t); % To access from places where handles.home_dir is unknown (must precede gateLoadFile())
if (~isempty(drv))
gateLoadFile(handles,drv,varargin{1}); % load recognized file types
else
recentFiles(handles,[]); % Just make the "Recent files" entry available
end
if (~isempty(handles.hImg)) % If we had something in input, check the coordinates type
handles = aux_funs('isProj', guidata(hObject));
setAxesDefCoordIn(handles,1);
end
try custom_menus(hObject, handles.path_data), end % Add any custom menus as commanded by OPTcontrol.txt
if (~isempty(feval_me)) % 'feval_me' must be a string with a command that calls a child figure.
try
c = false(numel(x_comm),1);
for (k = 1:numel(x_comm)) % Retain only eventual commands to be send to children
if (~isempty(x_comm{k})), c(k) = true; end
end
if (any(c)), x_comm = x_comm(c);
else, x_comm = {[]};
end
ind = strfind(feval_me, ',');
if (isempty(ind))
h = feval(feval_me, x_comm);
else
gui_Callback = str2func(feval_me(1:ind(1)-1));
arg = feval_me(ind(1)+1:end); % Restrict to one arg only because I have not yet use cases with more
ind2 = strfind(arg, '(');
if (isempty(ind2))
h = feval(gui_Callback, arg, x_comm); % What case is this ???
else % For example, guidata(gcf)
handles = guidata(handles.figure1); % Need an updated one
if (strcmp(arg, 'guidata(gcf)'))
h = feval(gui_Callback, handles, x_comm);
elseif (strcmp(arg, 'gcf'))
h = feval(gui_Callback, handles.figure1, x_comm);
else
if (strncmp(arg,'guidata(gcf),',13))
h = feval(gui_Callback, handles, arg(14:end), x_comm);
else
h = feval(gui_Callback, arg, x_comm); % Does this case exists ???
end
end
end
end
catch
disp(lasterr)
end
hObject = struct('hMirFig',hObject, 'hChildFig',h); % Because we may need the handles of both Figs
end
% --------------------------------------------------------------------------------------------------
function erro = gateLoadFile(handles,drv,fname)
% Gateway function to load a recognized file type using its name
erro = 0;
if (strncmp(drv, 'MB', 2)), drv = 'MB'; end % The MB (system) type case are actually several
switch drv
case 'gmt', loadGRID(handles, fname, 'GMT_relatives')
case 'dono', erro = FileOpenGeoTIFF_CB(handles,'dono',fname); % It means "I don't know" and dealt by GDAL
case 'generic', FileOpenNewImage_CB(handles, fname);
case 'geotif', FileOpenGeoTIFF_CB(handles, 'nikles', fname);
case 'ecw', FileOpenGeoTIFF_CB(handles, 'ecw', fname); % A particular case (includes jp2)
case 'multiband', FileOpenGDALmultiBand_CB(handles, 'AVHRR', fname);
case 'envherd', FileOpen_ENVI_Erdas_CB(handles, [], fname);
case 'mat', FileOpenSession_CB(handles, fname)
case 'mola', loadGRID(handles, fname, 'MOLA');
case 'cpt', color_palettes(fname);
case 'dat', load_xyz(handles, fname);
case 'pick', load_xyz(handles, fname, drv);
case 'ncshape', load_xyz(handles, fname, drv);
case 'ncchimoce', load_xyz(handles, fname, drv);
case 'shp', DrawImportShape_CB(handles, fname);
case 'ogr', DrawImportOGR_CB(handles, fname);
case 'las', read_las(handles, fname);
case 'mgg_gmt', GeophysicsImportGmtFile_CB(handles, fname);
case 'ghost', load_ps(handles, fname); % Put ghostscript on works
case 'sww', aquamoto(fname);
case 'MB', show_mb(handles, fname);
case 'xtf', showXTF_CB(handles, fname);
otherwise, erro = 1;
end
if (erro), warndlg(['Sorry but couldn''t figure out what to do with the ' fname ' file'],'Warning'), end
% --------------------------------------------------------------------------------------------------
function handles = recentFiles(handles, opt)
if (~isempty(handles.fileName) && ~exist(handles.fileName,'file')), return, end % For example, subdatasets
jump = false; N = min(numel(handles.FOpenList), numel(handles.RecentF));
if (nargin == 1 && ~isempty(handles.fileName)) % Update list
for (i = 1:N) % See if new name is already in FOpenList
if (strcmpi(handles.fileName, handles.FOpenList{i}))
for (k = i:-1:2) % Yes, it is. So move it to the top
handles.FOpenList{k} = handles.FOpenList{k-1};
end
handles.FOpenList{1} = handles.fileName;
jump = true; break
end
end
if (~jump) % We got a new name
for (i = N:-1:2) % Make room for the new name at top of the list
handles.FOpenList{i} = handles.FOpenList{i-1};
end
handles.FOpenList{1} = handles.fileName;
end
end
if (isempty(handles.fileName)), jump = true; end
if ( ~jump && (nargin == 1 || (nargin == 2 && ~isempty(opt))) ) % Save only if it worth it
FOpenList = handles.FOpenList; fname = [handles.path_data 'mirone_pref.mat'];
if (~handles.version7), save(fname,'FOpenList','-append') % Update the list for "Recent files"
else, save(fname,'FOpenList','-append', '-v6')
end
end
for (i = 1:N) % The ishandle test below is crutial when GCPs
if (isempty(handles.FOpenList{i}) || ~ishandle(handles.RecentF(i))), continue, end
set(handles.RecentF(i), 'Label',handles.FOpenList{i},'Callback',{@openRF,i}, 'Vis','on')
end
if (~nargout), guidata(handles.figure1,handles), end
% --------------------------------------------------------------------------------------------------
function openRF(obj,event,n)
% Open from 'Recent Files'
handles = guidata(obj);
[drv, sim] = aux_funs('findFileType',handles.FOpenList{n});
if (sim)
gateLoadFile(handles,drv,handles.FOpenList{n});
else % File does not exist. Update the FOpenList
warndlg(['File: ' handles.FOpenList{n} ' no longer exists'],'Warning')
handles.FOpenList(n) = []; handles.FOpenList{end+1} = []; % Delete missing and add one at the end
delete(handles.RecentF(n)), guidata(handles.figure1,handles)
end
% --------------------------------------------------------------------------------------------------
function handles = SetAxesNumericType(handles,event)
% Save original X & Y labels in appdata for easear access when we want to change them
setappdata(handles.axes1,'XTickOrig',get(handles.axes1,'XTickLabel'))
setappdata(handles.axes1,'XTickOrigNum',get(handles.axes1,'XTick'))
setappdata(handles.axes1,'YTickOrig',get(handles.axes1,'YTickLabel'))
setappdata(handles.axes1,'YTickOrigNum',get(handles.axes1,'YTick'))
n1 = numel(get(handles.axes1,'YTick'));
fsize = get(handles.axes1, 'FontSize');
set(handles.axes1, 'FontSize', 9) % Make this the default
n2 = numel(get(handles.axes1,'YTick'));
if (n1 ~= n2) % BUT SOMETIMES IT FCK IT??????????. DAMN ML BUGS
set(handles.axes1, 'FontSize', fsize)
end
LFT = 'DegDec'; visibility = 'on'; % For the geog case
if (~handles.geog), LFT = 'NotGeog'; visibility = 'off'; end
setappdata(handles.axes1,'LabelFormatType',LFT)
set(handles.LabFormat, 'Vis', visibility)
if (handles.validGrid), set(handles.PixMode, 'Callback', {@PixMode_CB,handles.figure1, true}, 'Vis', 'on')
else, set(handles.PixMode, 'Vis', 'off')
end
set(handles.RCMode, 'Callback', {@PixMode_CB,handles.figure1, false})
set(handles.FancyMode, 'Callback', {@FancyMode_CB,handles.figure1})
% --------------------------------------------------------------------------------------------------
function PixMode_CB(hObject, event, hFig, opt)
% Inside each grid cell, which is a pixel on the screen, display only the grid node value
handles = guidata(hFig);
if (opt) % Pixel mode on/off
if (strcmp(get(hObject,'Checked'),'off'))
set(hObject,'Checked','on'), setappdata(hFig,'PixelMode',true)
else
set(hObject,'Checked','off'), setappdata(hFig,'PixelMode',false)
end
set(handles.RCMode, 'Checked','off'), setappdata(hFig,'RCMode',false) % Put the RowColMode to off
else % Row Col mode on/off
if (strcmp(get(hObject,'Checked'),'off'))
set(hObject,'Checked','on'), setappdata(hFig,'RCMode',true)
else
set(hObject,'Checked','off'), setappdata(hFig,'RCMode',false)
end
set(handles.PixMode, 'Checked','off'), setappdata(hFig,'PixelMode',false) % Put the PixMode to off
end
% --------------------------------------------------------------------------------------------------
function FancyMode_CB(hObject, evt, hFig)
% Call the function that sets or unsets the FANCY frame depending on the previous state
handles = guidata(hFig);
if (strcmp(get(hObject,'Checked'),'off'))
set(hObject,'Checked','on'), fancyFrame(handles,'set')
else
set(hObject,'Checked','off'), fancyFrame(handles,'unset')
end
% --------------------------------------------------------------------------------------------------
function PlatesAgeLift_CB(handles)
% Apply Parsons & Sclatter relation to compensate sea-bottom age sinking.
% It assumes that bathymetry is the loaded grid (in meters Z up) and age in Ma
if (handles.no_file), return, end
resp = inputdlg({'Full name (with path) of Age grid:'},'Where is the grid?',[1 80]);
if (isempty(resp)), return, end
fname = resp{1};
if (~exist(fname,'file'))
warndlg('The file name provided does not exist. Bye, Bye','Warning'), return
end
att = gdalread(fname,'-M','-C');
if (att.GMT_hdr(1) > handles.head(1) || att.GMT_hdr(2) < handles.head(2) || ...
att.GMT_hdr(3) > handles.head(3) || att.GMT_hdr(4) < handles.head(4))
errordlg('No way. The Age grid does not cover the limits of the bathymetry grid.','Error'), return
end
[Age, X, Y, srsWKT, miniHandles] = read_grid([], fname, 'GDAL', sprintf('-R%.12f/%.12f/%.12f/%.12f', handles.head(1:4)));
if (isempty(Age)), return, end % Something bad happened
lift = single(350 * sqrt(double(Age)));
[X,Y,Z] = load_grd(handles);
lift = cvlib_mex('resize', lift, [size(Z,1) size(Z,2)]);
lifted = cvlib_mex('add', lift, Z);
miniHandles.head(7:9) = handles.head(7:9); % min/max (5:6) will be updated in GRDdisplay
GRDdisplay(handles,X,Y,lifted,miniHandles.head,'','AgeLiftedBathymetry',srsWKT)
% Now compute the depth anomaly
cvlib_mex('addS', lift, 2500);
cvlib_mex('add', lift, Z);
GRDdisplay(handles,X,Y,lift,miniHandles.head,'','Bathymetry anomaly',srsWKT)
% --------------------------------------------------------------------------------------------------
function load_ps(handles, fname)
% Load a PostScript file ... via ghostscript and a MEX to access ghost stdout (SHITY QUALITY)
% P6
% # Image generated by Ghostscript (device=ppmraw)
% 612 792
% 255
P = popenr(['gswin64c -q -r300x300 -sDEVICE=ppmraw -sOutputFile=- ' fname]);
try
% The following while's are because we can't fseek a pipe and header size is variable.
while (popenr(P,1,'char') ~= 10), end % Consume First header line
while (popenr(P,1,'char') ~= 10), end % Consume Second header line
t = ' ';
for (k = 1:20)
t(k) = popenr(P,1,'char');
if (t(k) == 10), break, end % EOF reached
end
t(k:end) = []; % Remove excess characters
ii = find(t == 32); % Find the dimensions separator (a blank)
width = str2double(char( t(1 : ii(1)-1) )');
height = str2double(char( t(ii(1)+1 : end) )');
while (popenr(P,1,'char') ~= 10), end % Consume Fourth header line
img = popenr(P, [height width 3], 'char');
catch
popenr(P, -1) % Close the stdin stream
errordlg(lasterr, 'Error'), return
end
popenr(P, -1) % Close the stdin stream
mirone(img);
% --------------------------------------------------------------------------------------------------
function load_MB(handles, fname, frmt)
% Send MB datalist or single MB file to mbimport and get back an Image like that produced by mbswath
cmd = '';
fid = fopen(fname,'rt');
todos = fread(fid,'*char'); fclose(fid);
if (todos(1) == '#' || todos(1) == '>')
ind = find(todos == 10); % Find new lines
cmd = todos(2:ind(1)-1)';
%todos = todos(ind(1)+1:end);
end
%[nomes, frmt] = strread(todos,'%s %s');
I = gmtmex(['mbimport -I' fname ' ' cmd]);
mirone(I)
if (handles.no_file), delete(handles.figure1), end % Delete the old and empty fig.
% --------------------------------------------------------------------------------------------------
function varargout = ImageCrop_CB(handles, opt, opt2, opt3)
% OPT is either a handle to a line that may be a rectangle/polygon,
% a Mx2, a 2xM matrix with the (x y) vector coordinates, or a [xmin xmax ymin ymax] BB.
% If empty, it calls rubberbandbox to do a interactive croping (called by "Crop Grid")
% OPT2 is a string to direct this function to different operations that
% apply to the grid and update the image.
% OPT3 contains the interpolation method when OPT2 == 'FillGaps' ('cubic', 'linear' or 'sea')
% Or the fill value when OPT2 == 'SetConst'
% Note: I won't make the "Drape" option active in the cropped window
%
% VARARGOUT -> If used will hold the result of this function instead of creating a new Fig
% Currently implemented in cases:
% Crop image (opt == hLine), 'CropaWithCoords', 'CropaGrid_pure', 'ROI_Mean', 'ROI_Median', 'ROI_STD'
% For the 'CropaGrid_pure' case varargout actualy returns 4 values -- {X,Y,Z_rect,head}
if (handles.no_file), return, end
first_nans = false; pal = []; mask = []; crop_pol = false; % Defaults to croping from a rectangle
wasROI = false; done = false; invert = false;
if (nargin < 3), opt2 = []; end
if (nargin < 4), opt3 = []; end
if ~isempty(opt) % OPT must be a rectangle/polygon handle (the rect may serve many purposes)
if ((numel(opt) == 1) && ishandle(opt))
x = get(opt,'XData'); y = get(opt,'YData');
elseif (numel(opt) == 4) % Assume that we have a [xmin xmax ymin ymax] BB
x = [opt(1) opt(1) opt(2) opt(2) opt(1)];
y = [opt(3) opt(4) opt(4) opt(3) opt(3)];
else
if (size(opt,2) > 2), x = opt(1,1:end); y = opt(2,1:end); % Row vectors
else, x = opt(:,1)'; y = opt(:,2)'; % Were col vectors, make them row for consistency
end
end
xp = [min(x) max(x)];
yp = [min(y) max(y)];
if (numel(x) == 5 && (x(1) == x(end)) && (y(1) == y(end)) && ... % Test if we have a rectangle
(x(1) == x(2)) && (x(3) == x(4)) && (y(1) == y(4)) && (y(2) == y(3)) )
if (strcmp(opt2,'SetConst'))
[xp, yp] = aux_funs('adjust_rect', handles, xp, yp); % Adjust such that only inside nodes are selected
end
else
dx2 = ~handles.head(7) * handles.head(8)/2; dy2 = ~handles.head(7) * handles.head(9)/2;
if (xp(1) < handles.head(1)-dx2 || xp(2) > handles.head(2)+dx2 || yp(1) < handles.head(3)-dy2 || yp(2) > handles.head(4)+dy2)
% Somewhat rare case where the polygon extends to outside grid/img limits. Must crop it to them.
P1.x = x(:); P1.y = y(:); P1.hole = 0; P2.hole = 0;
P2.x = [handles.head(1)-dx2; handles.head(1)-dx2; handles.head(2)+dx2; handles.head(2)+dx2; handles.head(1)-dx2];
P2.y = [handles.head(3)-dy2; handles.head(4)+dy2; handles.head(4)+dy2; handles.head(3)-dy2; handles.head(3)-dy2];
outPolyg = PolygonClip(P1, P2); % Intersection of polygon and map limits
x = outPolyg.x; y = outPolyg.y;
xp(1) = min(x); xp(2) = max(x);
yp(1) = min(y); yp(2) = max(y);
end
x_lim = [min(x) max(x)]; y_lim = [min(y) max(y)];
crop_pol = true; % Flag that we are croping from a polygon
end
rect_crop = [xp(1) yp(1) (xp(2) - xp(1)) (yp(2) - yp(1))]; % Even when from a rectangle (to insure a known order)
if (isempty(opt2) || strcmp(opt2,'CropaWithCoords')) % Pure Image croping OR Crop Image with coordinates
Z_rect = get(handles.hImg,'CData');
if (isa(Z_rect, 'uint16')), Z_rect = scaleto8(Z_rect); end % In the rare event of "Noe Diluge" images
[I,r_c] = cropimg(handles.head(1:2),handles.head(3:4),Z_rect,rect_crop,'out_grid');
if (crop_pol) % Shape cropping
x_lim(1) = handles.head(1) + (r_c(3)-1)*handles.head(8); x_lim(2) = handles.head(1) + (r_c(4)-1)*handles.head(8);
y_lim(1) = handles.head(3) + (r_c(1)-1)*handles.head(9); y_lim(2) = handles.head(3) + (r_c(2)-1)*handles.head(9);
mask = ~(img_fun('roipoly_j',x_lim,y_lim,I,x,y));
if (ndims(I) == 2)
I(mask) = 0;
else
for (k = 1:3)
tmp = I(:,:,k); tmp(mask) = 255; I(:,:,k) = tmp;
end
mask = uint8(~mask); cvlib_mex('CvtScale',mask,255,0); % NEW. 1st transparency attempt
I(:,:,4) = mask;
clear tmp mask
end
end
[m,n] = size(I);
elseif (strcmp(opt2,'ImplantGrid')) % Read and external grid and implant it in host grid
[Z_rect, r_c] = transplants(opt, 'grid', true, handles);
if (isempty(Z_rect)), return, end % User gave up
% elseif (strcmp(opt2,'CropaStack')) % Crop a multi-layered file
% BL = getappdata(handles.figure1,'BandList');
% n_layers = BL{7}(end);
% Z_rect = aux_funs('get_layer_n', handles.figure1, 1);
% [I,r_c] = cropimg(handles.head(1:2),handles.head(3:4),Z_rect,rect_crop,'out_grid');
% I = repmat(I, [1, 1, n_layers]); % Make room for the remaining layers
% for (k = 2:n_layers)
% Z_rect = aux_funs('get_layer_n', handles.figure1, k);
% [I(:,:,k),r_c] = cropimg(handles.head(1:2),handles.head(3:4),Z_rect,rect_crop,'out_grid');
% end
else % Extract the sub-grid inside the rectangle/polygon
[X,Y,Z,head] = load_grd(handles);
if isempty(Z), return, end % An error message was already issued
[Z_rect,r_c] = cropimg(head(1:2),head(3:4),Z,rect_crop,'out_grid');
if (crop_pol)
zzz = grdutils(Z_rect,'-L'); z_min = zzz(1); clear zzz;
resp = [];
if (strncmp(opt2,'CropaGrid_pure', 14))
if (opt2(end) == '-') % if we get a 'CropaGrid_pure-' don't ask for the fill value. Just use NaNs
resp = NaN;
opt2(15:end) = []; % Remove so not to perturb later tests
else
resp = inputdlg({'Enter outside polygon value'},'Choose out value',[1 30],{sprintf('%.4f',z_min)}); pause(0.01)
if isempty(resp), return, end
resp = str2double(resp{1});
end
elseif (strcmp(opt2,'ROI_SetConst')) % Set the polygon in-or-out to cte
resp = question({'Enter new grid value'},'Replace with cte value',[1 30],'NaN','whatever');
if isempty(resp), return, end
invert = resp{2};
resp = str2double(resp{1});
end
if (isnan(resp)), handles.have_nans = 1; end
x_lim(1) = head(1) + (r_c(3)-1)*head(8); x_lim(2) = head(1) + (r_c(4)-1)*head(8);
y_lim(1) = head(3) + (r_c(1)-1)*head(9); y_lim(2) = head(3) + (r_c(2)-1)*head(9);
mask = img_fun('roipoly_j', x_lim, y_lim, Z_rect, x, y);
if (strcmp(opt2,'CropaGrid_pure'))
Z_rect(~mask) = single(resp);
elseif (strcmp(opt2,'ROI_SetConst'))
if (invert), mask = ~mask; end % Mask in the outside instead
Z_rect(mask) = single(resp); % Set the mask values to const
if (invert)
Z(:) = single(resp); % Actually, we need to mask the whole outside
else
handles.Z_back = Z(r_c(1):r_c(2),r_c(3):r_c(4)); handles.r_c = r_c; % For the Undo op
end
Z(r_c(1):r_c(2),r_c(3):r_c(4)) = Z_rect;
if (isnan(resp)), handles.have_nans = 1; first_nans = true; end
elseif (strcmp(opt2,'ROI_Mean'))
res = Z_rect(mask); res(isnan(res)) = []; res = mean(double(res));
if (nargout), varargout{1} = res; return, end
msgbox(sprintf('Mean over ROI = %15g',res),'ROI-Mean'), return
elseif (strcmp(opt2,'ROI_Median'))
res = Z_rect(mask); res(isnan(res)) = []; res = median(res);
if (nargout), varargout{1} = res; return, end
msgbox(sprintf('Median over ROI = %15g',res),'ROI-Median'), return
elseif (strcmp(opt2,'ROI_STD'))
res = Z_rect(mask); res(isnan(res)) = []; res = std(double(res));
if (nargout), varargout{1} = res; return, end
msgbox(sprintf('STD over ROI = %15g',res),'ROI-STD'), return
elseif (strcmp(opt2,'ROI_MedianFilter'))
[Z,Z_rect,handles] = roi_filtering(handles, Z, head, Z_rect, r_c, mask);
elseif (strcmp(opt2,'ROI_SplineSmooth'))
opt2 = 'SplineSmooth'; % Strip the 'ROI_' part so that we can use the same code as for rectangles
wasROI = true; % Signal the SplineSmooth code below that we need to mask result
elseif (strncmp(opt2,'ROI_FillSinks', 13))
opt2 = opt2(5:end); % Strip the 'ROI_' part so that we can use the same code as for rectangles
wasROI = true; % Signal the FillSinks_xxxx code below that we need to mask result
elseif (strcmp(opt2,'CropaGrid_histo'))
Z_rect(~mask) = single(NaN);
elseif (strcmp(opt2,'ROI_Clip'))
opt2 = 'Clip'; % Now that we have the mask, make this case == to rectangle clip
else
errordlg('Unknown case in ImageCrop','Error'), return
end
end
[m,n] = size(Z_rect);
end
else % Interactive croping (either Grid or Image)
if (strcmp(opt2,'CropaGrid')) % Arrive here when called by "Grid Tools -> Crop Grid"
[X,Y,Z,head] = load_grd(handles);
if isempty(Z), return, end
[p1,p2] = rubberbandbox;
x0 = min(p1(1),p2(1)); y0 = min(p1(2),p2(2));
dx = abs(p2(1)-p1(1)); dy = abs(p2(2)-p1(2));
[Z_rect,r_c] = cropimg([head(1) head(2)],[head(3) head(4)],Z,[x0 y0 dx dy],'out_grid');
X = linspace( head(1) + (r_c(3)-1)*head(8), head(1) + (r_c(4)-1)*head(8), r_c(4) - r_c(3) + 1 );
Y = linspace( head(3) + (r_c(1)-1)*head(9), head(3) + (r_c(2)-1)*head(9), r_c(2) - r_c(1) + 1 );
head(1) = X(1); head(2) = X(end); head(3) = Y(1); head(4) = Y(end);
GRDdisplay(handles,X,Y,Z_rect, head, '', sprintf('Cropped_(%.0f)_grid', rand(1)*100))
return
else % Just a image crop op
I = cropimg; [m,n] = size(I);
end
end
if (isempty(opt2) || strcmp(opt2,'CropaWithCoords')) % Just pure Image croping
if (m < 2 || n < 2), return, end % Image too small.
if (strcmp(get(handles.axes1,'Ydir'),'normal')), I = flipdim(I,1); end