; Default basic data info. codeName = 'cai_map-2sep16_DE.pro of 2-Sep-2016 6pm' TRUE = 1 FALSE = 0 maxWhite = 255 ; white value in tiff images. valBlk = 0 beamCurrent = 0.0 tempul = ULONG(0) ; temporary value to use as needed. tempul2 = ULONG(0) ; nColors = 17 ; nPhases = 17 ; Declare a set of indices for the array of elements indxXx = 0 ;the Mask indxMg = 1 indxCa = 2 indxAl = 3 indxSi = 4 indxTi = 5 indxFe = 6 indxS = 7 indxNi = 8 indxNa = 9 indxBE = 10 sizeElmArray = 11 ; Declare a 'structure' with the name 'elm', with element information. ; items are: index value, name, atom # z, present? T/F, index in MCA maps, index in STF maps, count(intensity sum) struct = { elm, index:0, present:0, label:'', z:0,iMCA:0, iSTF:0, count:0L } ; Make an array of structures with the values of 'elm' given by the values in the array. elms = REPLICATE({elm}, sizeElmArray) elms[indxXx ] = {elm, 0, 1, 'mask', 0, 0, 0, 0L } elms[indxMg ] = {elm, 1, 1, 'Mg', 12, 1, 0, 0L } elms[indxCa ] = {elm, 2, 1, 'Ca', 20, 2, 0, 0L } elms[indxAl ] = {elm, 3, 1, 'Al', 13, 3, 0, 0L } elms[indxSi ] = {elm, 4, 1, 'Si', 14, 0, 1, 0L } elms[indxTi ] = {elm, 5, 1, 'Ti', 22, 0, 2, 0L } elms[indxFe ] = {elm, 6, 1, 'Fe', 26, 0, 3, 0L } elms[indxS ] = {elm, 7, 1, 'S', 16, 0, 0, 0L } elms[indxNi ] = {elm, 8, 1, 'Ni', 28, 0, 0, 0L } elms[indxNa ] = {elm, 9, 1, 'Na', 11, 0, 0, 0L } elms[indxBE ] = {elm,10, 0, 'BSE', 1, 0, 0, 0L } ; Declare a set of easy-to-read index values that will be used to index the values in the array 'cat[n]' ; The mask background in result will be white, the unknown interior light gray (ltgray). indxMsk = 0 ; mask indxUnk = 1 ; unknown indxHole = 2 ; nothing indxCc = 3 ; carbonate indxMet = 4 ; metal indxFeS = 5 ; troilite-FeS indxCaS = 6 ; CaS indxMnS = 7 ; MnS indxMgS = 8 ; MgS indxSpn = 9 ; MgAl2O4-spinel indxTO2 =10 ; TiO2 indxPrv =11 ; perovskite-CaTiO3 indxTtn =12 ; sphene-CaTiSiO5 indxSiO2 =13 ; silica indxWol =14 ; wollastonite-CaSiO3 indxCor =15 ; corundum indxHib =16 ; hibonite indxCA2 =17 ; grossite-CA2 indxCA1 =18 ; CaAl2O4-CA1 indxOlv =19 ; olivine indxGls =20 ; glass indxOpx =21 ; Ca-poor-pyroxene indxCpx =22 ; Ca-rich-pyroxene indxMel =23 ; melilite indxPlg =24 ; anorthite indxAlt =25 ; Na-silicate-alteration indxRhn =26 ; Rhonite indxIlm =27 ; ilmenite-FeTiO3 indxAb =28 ; Albite indxHrc =29 ; Hercynite indxCaTs =30 ; CaAlAlSiO6 indxMsc =31 ; miscellaneous sizePhaseArray = 32 ; Declare the indices for certain ratios: indxRatCaAl = 0 indxRatTiSi = 1 indxRatCaSi = 2 indxRatMgSi = 3 indxRatFeSi = 4 ; Declare a 'structure' with the name 'phase', consisting of various items also declared here. ; so 'count' will be a long integer, initialized with a value of 0. And 'rgbX' is an integer array for an RGB color. ; 25jan2012: Add rgbMCALim1 and rgbSTFLim2: These are max and min RGB values for selection criteria in Mg-Ca-Al map, and then for secondary map. ; 19feb2012: Make 2nd set of 3 for Si-Ti-Fe limits. ; 2mar2012: Add elemMeans arrays. ; 12jul2016: Make elemSums a DBLARR (DSE) struct = { phase, index:0, gsVal1:0, label:'', count:0L, color:'', rgbX:INTARR(3), phsIncl:0, gsVal2:0, $ ratCrit:FLTARR(5), ratRange:FLTARR(5), elemSums:DBLARR(10), elemMeans:FLTARR(10)} ; gray = 204,204,204 ;--------------------------i--gsVal1--------- phase name ------------ RGB color ---- RGB colors --include-gs2--- ratio Criteria -------- ratio Ranges ------------------- element sums ------------------------- element mean values --------- ; Make an array of structures with the values of 'phase' given by the values in the array. ; 1 2 3 4 5 --- 1 2 3 4 5 --- 1 2 3 4 5 6 7 8 9 10 --- 1 2 3 4 5 6 7 8 9 10 cat = REPLICATE({phase}, sizePhaseArray) cat[indxMsk ] = {phase, 0, 255 , 'mask', 0, 'RGBwhite', [ 255, 255, 255], TRUE, 255, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxUnk ] = {phase, 1, 247 , 'unknown', 0, 'RGBblack', [ 0, 0, 0], TRUE, 200, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxHole ] = {phase, 2, 254 , 'hole', 0, 'white', [ 254, 254, 254], TRUE, 0 , [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxCc ] = {phase, 3, 229 , 'carbonate', 0, 'brown', [ 96, 57, 19], TRUE, 40 , [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxMet ] = {phase, 4, 179 , 'metal', 0, 'dark-grey', [ 65, 64, 66], TRUE, 90 , [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxFeS ] = {phase, 5, 165 , 'troilite-FeS', 0, 'lt-brown', [ 169, 124, 80], TRUE, 100, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxCaS ] = {phase, 6, 160 , 'CaS', 0, 'lt-yellow', [ 246, 236, 149], TRUE, 240, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxMnS ] = {phase, 7, 162 , 'MnS', 0, 'lt-yellow', [ 246, 236, 149], TRUE, 242, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxMgS ] = {phase, 8, 167 , 'MgS', 0, 'yellow', [ 247, 229, 48], TRUE, 244, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxSpn ] = {phase, 9, 25 , 'MgAl2O4-spinel', 0, 'purple', [ 127, 71, 155], TRUE, 130, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxTO2 ] = {phase, 10, 60 , 'TiO2', 0, 'grey', [ 147, 149, 152], TRUE, 246, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxPrv ] = {phase, 11, 125 , 'perovskite-CaTiO3', 0, 'brit-yellow', [ 253, 215, 0], TRUE, 145, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxTtn ] = {phase, 12, 60 , 'sphene-CaTiSiO5', 0, 'pale-blue', [ 162, 216, 232], TRUE, 248, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxSiO2 ] = {phase, 13, 125 , 'silica', 0, 'lavender', [ 152, 153, 203], TRUE, 15, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxWol ] = {phase, 14, 60 , 'wollastonite-CaSiO3', 0, 'orange', [ 232, 116, 38], TRUE, 250, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxCor ] = {phase, 15, 60 , 'corundum', 0, 'light-green', [ 177, 197, 54], TRUE, 245, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxHib ] = {phase, 16, 115 , 'hibonite', 0, 'hot-pink', [ 230, 65, 151], TRUE, 170, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxCA2 ] = {phase, 17, 102 , 'grossite-CA2', 0, 'magenta', [ 158, 31, 99], TRUE, 180, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxCA1 ] = {phase, 18, 76 , 'CaAl2O4-CA1', 0, 'pink', [ 242, 119, 152], TRUE, 210, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxOlv ] = {phase, 19, 217 , 'olivine', 0, 'red', [ 237, 28, 36], TRUE, 50 , [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxGls ] = {phase, 20, 204 , 'glass', 0, 'navy-blue', [ 54, 55, 149], TRUE, 60 , [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxOpx ] = {phase, 21, 191 , 'Ca-poor-pyroxene', 0, 'dark-red', [ 151, 51, 51], TRUE, 70 , [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxCpx ] = {phase, 22, 50 , 'Ca-rich-pyroxene', 0, 'green', [ 44, 126, 62], TRUE, 120, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxMel ] = {phase, 23, 153 , 'melilite', 0, 'lite-orange', [ 248, 152, 57], TRUE, 150, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxPlg ] = {phase, 24, 88 , 'anorthite', 0, 'steel-blue', [ 28, 117, 188], TRUE, 190, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxAlt ] = {phase, 25, 140 , 'Na-alteration', 0, 'dark-blue', [ 33, 64, 154], TRUE, 125, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxRhn ] = {phase, 26, 185 , 'Rhonite', 0, 'sea-green', [ 57, 181, 74], TRUE, 135, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxIlm ] = {phase, 27, 60 , 'ilmenite-FeTiO3', 0, 'light-grey', [ 209, 211, 212], TRUE, 248, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxAb ] = {phase, 28, 85 , 'albite', 0, 'teal', [ 36, 184, 148], TRUE, 192, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxHrc ] = {phase, 29, 25 , 'hercynite', 0, 'light-pink', [ 246, 151, 153], TRUE, 140, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxCaTs ] = {phase, 30, 50 , 'Ca-Alpyroxene', 0, 'aqua', [ 9, 194, 236], TRUE, 115, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } cat[indxMsc ] = {phase, 31, 12 , 'miscellaneous', 0,'forest-green', [ 24, 63, 29], TRUE, 260, [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0], [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0] } ; NOTES: in older CAImode3, there is a different set of 'modeRGBcol' which are values in 'rgbX' in 'cat' array. ; Set a few more variables that will be useful ; OUTPUT mask will be white in final output. INPUT mask is BLACK.!! countSul = ULONG(0) ; count total of all sulfide phases. gsMsk = cat[indxMsk].gsVal1 ; grayscale value of the mask. (WHITE) gsUnk = cat[indxUnk].gsVal1 ; grayscale value of the unknowns. ; Set up counters for pixels total, and object. countObj = 0L countTot = 0L elemAvgVal = DBLARR(sizeElmArray, 6) ; change to DBL 12-July-2016 DSE (max FLOAT is ~2.1 billion) Add a value also. ;*************************************************************************************** ; Constants for application of Bence-Albee correction (Clarke et al. 2001, J Pet) ; alpha factors: Love & Scott alpha factors, from Armstrong 1984 'Microbeam Analysis' ; Si02 Al203 FeO MnO MgO CaO Na20 K20 TiO2 NiO SO2 Sc2O3 V2O3 ; Si 1.000 1.494 1.279 1.211 1.431 1.048 1.356 0.987 1.000 1.000 1.000 ; Al 0.999 1.000 1.556 1.451 1.756 1.149 1.645 1.072 1.000 1.000 1.000 ; Fe 1.148 1.128 1.000 0.987 1.141 1.124 1.114 1.093 1.000 1.000 1.000 ; Mn 1.161 1.140 1.014 1.000 1.153 1.149 1.125 1.117 1.000 1.000 1.000 ; Mg 1.086 1.009 2.013 1.846 1.000 1.302 2.134 1.196 1.000 1.000 1.000 ; Ca 1.091 1.068 0.932 0.902 1.073 1.000 1.044 1.197 1.000 1.000 1.000 ; Na 1.352 1.242 3.066 2.778 1.130 1.741 1.000 1.564 1.000 1.000 1.000 ; K 1.144 1.118 0.987 0.959 1.120 0.858 1.088 1.000 1.000 1.000 1.000 ; Ti 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ; Ni 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ; S 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ; Cr 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ; Sc 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 ; V 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 alpha = FLTARR(14,14) ; Si Al Fe Mn Mg Ca Na K Ti Ni S Cr Sc V alpha = [ [1.000,1.494,1.279,1.211,1.431,1.048,1.356,0.987,1.000,1.000,1.000,1.000,1.000,1.000],$ [0.999,1.000,1.556,1.451,1.756,1.149,1.645,1.072,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.148,1.128,1.000,0.987,1.141,1.124,1.114,1.093,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.161,1.140,1.014,1.000,1.153,1.149,1.125,1.117,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.086,1.009,2.013,1.846,1.000,1.302,2.134,1.196,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.091,1.068,0.932,0.902,1.073,1.000,1.044,1.197,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.352,1.242,3.066,2.778,1.130,1.741,1.000,1.564,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.144,1.118,0.987,0.959,1.120,0.858,1.088,1.000,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000],$ [1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000,1.000] ] ; Note IDL uses column-row notation (like images, NOT mathematics). ; So a[3,2] refers to item in column 4, row 3. ; beta factors: will be recalculated from alpha factors & std. oxide data beta = FLTARR(11) ; standard counts (From ECP_Stds_10july11a_0asciifile.txt, in speccaliinfo10July.xls ; standard info: Counts - peak-background in counts/second (K, Ni, S are dummy values of 1000) ; Si Al Fe Mn Mg Ca Na K Ti Ni S Cr Sc V stdcps = [ 8812, 18938, 5232, 2938, 11924, 4850, 1504, 1000, 18218, 1000, 1000, 5940, 16480, 20018 ] ; standard info: Oxide concentrations in wt% oxide: ; Si Al Fe Mn Mg Ca Na K Ti Ni S Cr Sc V ; WakefieldDI Al2O3 FayalB6 Rhodonit 273Forst CaSiO3 albite K-min CaTiO3 NiS FeS MgCr2O4 Sc CsV2O5 stdoxidefrac = [ 55.49, 100.0, 29.49, 45.23, 57.29, 48.28, 11.82, 10.0, 48.76, 10.0, 10.0, 79.04, 100.0, 10.0 ] ;*************************************************************************************** ; Flags: Set default values: boxScoreyes1no2 = 2 ; For std grids, a box inside the larger image will be averaged. delBrightPxlyes1no2 = 1 ; Set ultra-bright pixels to average neibhbor values. brightDiffCrit = 100.0 brightThreshold = 2000.0 fileSpecEnd = '' pxlIncrement1 = 1 ; To speed up tests, change from 1 to a larger number. Will default to 1 if 0. ECPmac1orDSEpc2 = 1 ; switch for file locations etc. (DSE 26feb) inputFromListYes1no2 = 1 ; if yes, then jumps to input segment for batch processing of a lot of maps. bypassModeMappingYes1no2 = 2 ; if yes, NO modes will be calculated (faster!) writeIntermedTifYes1no2 = 1 ; if no, don't print a lot of the tif output. flagMaskDifferYes1no2 = 1 ; Add 19-July-2016 DSE cropMapsyes1no2 = 1 findForstYes1no2 = 2 ; add 4-Sep-2016 DSE to find forsterite pixels. ;********************************************************************************************** ; Start of main setting of file location and input parameters - For treating single maps - if (inputFromListYes1no2 EQ 2) then begin ; Set call to the data on each probe map. mapSet = 17 CASE mapSet OF 17: BEGIN ; 19-July-2016 To run masked standard mosaic: One set of maps, 7 masks! if (ECPmac1orDSEpc2 EQ 1) then begin endif else begin sourceDrive = 'USB30FD' sourceBaseDir = sourceDrive + '\trial1-Oct16\' targetBaseDir = sourceDrive + '\trial1-Oct16\' ;logfile = targetBaseDir + 'logFile1.txt' logfile1 = targetBaseDir + 'batchlogphaseFile.txt' logfile2 = targetBaseDir + 'batchlogelementFile.txt' logfile3 = targetBaseDir + 'batchlogphseleFile.txt' logfile4 = targetBaseDir + 'batchlogphseleFile.txt' endelse flagMaskDifferYes1no2 = 1 ; Add 19-July-2016 DSE basefile = 'Colony-4595-1-t3-c1_' targetfile = 'Colony-4595-1-t3-c1_3Sept_' title = 'Colony-4595-1-t3-c1_3Sept_' fileSuffix = '' micronPerPxl = 1.0 ; micron per pixel beamCurrent = 50.0 normBeamCurrent = TRUE initMaskCol = 0.0 dwellTime = 30.0 ; milliseconds (ms) (all 20nA, r1-2 = 20ms, 1micron, r3 = 20ms, 2mu, r4 = 12ms, 2mu) normDwellTime = TRUE fileSpecs = ['xx', 'Mg', 'Al', 'Ca', 'Si', 'Ti', 'Fe'] ; Fixed order for elements & maps: >> 'xx', 'Mg', 'Ca', 'Al', 'Si', 'Ti', 'Fe', (not present: 'S', Ni') numImages = 0 fileSpecEnd = '' FOR i= 0, (sizeElmArray-1) DO elms[i].present = 0 FOR i= 0, indxS DO BEGIN elms[i].present = 1 numImages = numImages + 1 ENDFOR ; if numImages NE 9 then print, 'ERROR! Wrong # of elements !! ' noMask = FALSE ; This value (inputMaskGS) is necessary if the mask was not set up correctly. ; inputMaskGS = 5 MapCrit = 2 ; for 32-bit 60ms/pxl 20nA data. inBitDepth = 32 boxScoreyes1no2 = 2 delBrightPxlyes1no2 = 1 brightDiffCrit = 100.0 brightThreshold = 2000.0 pxlIncrement1 = 1 ; To speed up tests, change from 1 to a larger number. Will default to 1. cropMapsyes1no2 = 2 ; To crop the maps all to the mask dimensions plus 'cropBufferMapPxls'. cropBufferMapPxls = 0 ; Buffer = 3 pixels all around. END ELSE: print, 'mapSet has illegal value.' ENDCASE ; Open log file for single map: openw, loglun1, logfile1, error=errlog, /get_lun openw, loglun2, logfile2, error=errlog, /get_lun openw, loglun3, logfile3, error=errlog, /get_lun IF (errlog NE 0) then PRINT, 'ERRROR opening logfile', !ERROR_STATE.MSG printf, loglun1, 'code: ', STRING(9B), codeName printf, loglun2, 'code: ', STRING(9B), codeName printf, loglun3, 'code: ', STRING(9B), codeName printf, loglun1, 'date: ', STRING(9B), SYSTIME() printf, loglun2, 'date: ', STRING(9B), SYSTIME() printf, loglun3, 'date: ', STRING(9B), SYSTIME() printf, loglun1, '=====================================' printf, loglun2, '=====================================' printf, loglun3, '=====================================' ; BATCH INPUT *************************************************************************************************************** endif else begin ; START of segment for batch input 'inputFromListYes1no2' = 1 ; Default switches: inBitDepth = 32 cropMapsyes1no2 = 2 boxScoreyes1no2 = 2 delBrightPxlyes1no2 = 1 brightDiffCrit = 100.0 brightThreshold = 2000.0 pxlIncrement1 = 1 ; To speed up tests, change from 1 to a larger number. Will default to 1. cropMapsyes1no2 = 1 ; To crop the maps all to the mask dimensions plus 'cropBufferMapPxls'. cropBufferMapPxls = 5 ; Buffer = N pixels all around. ; goto, bypassECPfileInput ;*********************************************************************************************ECP segment ; THIS is for ECP data: flagMaskDifferYes1no2 = 2 ; Add 19-July-2016 DSE ;set phases on or off for p=1, (sizePhaseArray-1) do cat[p].phsIncl = FALSE ; ECP phase list: Spn/TO2/Prv/Ttn/SiO2/Wol/Cor/Hib/Olv/Opx/Cpx/Mel/Plg/Ilm/Ab cat[indxMsk ].phsIncl = TRUE cat[indxUnk ].phsIncl = TRUE cat[indxHole ].phsIncl = TRUE ; Set holes = mask color-1 (change from mask color 12-July-2016 DSE). cat[indxCc ].phsIncl = FALSE cat[indxMet ].phsIncl = TRUE cat[indxFeS ].phsIncl = TRUE cat[indxCaS ].phsIncl = FALSE cat[indxMnS ].phsIncl = FALSE cat[indxMgS ].phsIncl = FALSE cat[indxSpn ].phsIncl = TRUE cat[indxTO2 ].phsIncl = TRUE ;TiO2 cat[indxPrv ].phsIncl = TRUE ;CaTiO3 cat[indxTtn ].phsIncl = FALSE ;sphene-CaTiSiO5 = titanite cat[indxSiO2 ].phsIncl = TRUE cat[indxWol ].phsIncl = FALSE cat[indxCor ].phsIncl = TRUE cat[indxHib ].phsIncl = TRUE cat[indxCA2 ].phsIncl = FALSE cat[indxCA1 ].phsIncl = FALSE cat[indxOlv ].phsIncl = TRUE cat[indxGls ].phsIncl = TRUE ; mesostasis cat[indxOpx ].phsIncl = TRUE cat[indxCpx ].phsIncl = TRUE cat[indxMel ].phsIncl = TRUE cat[indxPlg ].phsIncl = TRUE cat[indxAlt ].phsIncl = FALSE cat[indxRhn ].phsIncl = FALSE cat[indxIlm ].phsIncl = TRUE cat[indxAb ].phsIncl = TRUE cat[indxHrc ].phsIncl = FALSE cat[indxCaTs ].phsIncl = FALSE cat[indxMsc ].phsIncl = FALSE ; inputFromListYes1no2 = YES ; SEE F:\mets-111221\00-TomoTools ... tomo_bat_proc2.pro ; FILE locations (change 20-Jun-2016 from format of 28-Aug-2012): if (ECPmac1orDSEpc2 EQ 1) then begin sourceDrive = '' csvInputFile = sourceDrive + '/Volumes/USB30FD/kainsaz22/kainsaz22.csv' sourceBaseDir = sourceDrive + '/Volumes/USB30FD/mode input/' maskFileDir = sourceDrive + '/Volumes/USB30FD/mask files/' sourceELMmapSubDir = sourceDrive + '/Volumes/USB30FD/mode input/' sourceSubDir = '' logfile1 = '/Volumes/USB30FD/kainsaz22/batchlogphaseFile.txt' logfile2 = '/Volumes/USB30FD/kainsaz22/batchlogelementFile.txt' logfile3 = '/Volumes/USB30FD/kainsaz22/batchlogphseleFile.txt' logfile4 = '/Volumes/USB30FD/kainsaz22/batchlogprogressFile.txt' targetBaseDir = sourceDrive + '/Volumes/USB30FD/kainsaz22/' maskedFileDir = sourceDrive + '/Volumes/USB30FD/kainsaz22/Masked Output/' CritMapDir = sourceDrive + '/Volumes/USB30FD/kainsaz22/Other Output/' ModeDir = sourceDrive + '/Volumes/USB30FD/kainsaz22/Mode/' endif else begin inputFromListTask = 1 ; task 1 = get element map I/pxl for LA-ICP masked areas. sourceDrive = 'E' task = 1 CASE task OF 1: BEGIN ; For ECP mode of whole objects file work sourceBaseDir = sourceDrive + ':\ECP-data\' targetBaseDir = sourceBaseDir csvInputFile = sourceBaseDir + '\TestBatchMode-13sep16-chons.csv' logfile1 = targetBaseDir + 'ECP_batchMode-15sep16-1phases.txt' logfile2 = targetBaseDir + 'ECP_batchMode-13sep16-2elsTotal.txt' logfile3 = targetBaseDir + 'ECP_batchMode-13sep16-3elsByPhase.txt' logfile4 = targetBaseDir + 'ECP_batchMode-13sep16-4logprogress.txt' maskFileDir = sourceBaseDir + 'mask_files\' sourceELMmapSubDir = sourceBaseDir + 'mode_input\' CritMapDir = targetBaseDir + '13sep16-newC\' ModeDir = targetBaseDir + '13sep16-newC\' maskedFileDir = targetBaseDir + '13sep16-newC\' END ELSE: print, 'TASK has illegal value.' ENDCASE endelse inBitDepth = 32 micronPerPxl = 1.0 ; micron per pixel criterionBeamCurrent = 20.0 ; standard nA for normalization normBeamCurrent = TRUE criterionDwellTime = 60.0 normDwellTime = TRUE initMaskCol = 0.0 ; mask is BLACK delBrightPxlyes1no2 = 1 brightDiffCrit = 100.0 brightThreshold = 2000.0 cropMapsyes1no2 = 1 bypassModeMappingYes1no2 = 2 writeIntermedTifYes1no2 = 1 ; if no, don't print a lot of the tif output. flagMaskDifferYes1no2 = 1 ; Add 19-July-2016 DSE bypassECPfileInput: ; Open log files for batch of maps: openw, loglun1, logfile1, error=errlog, /get_lun ; logfile1 will be OUTPUT TABLE for MODE openw, loglun2, logfile2, error=errlog, /get_lun ; logfile2 will be OUTPUT TABLE for ELEMENTS openw, loglun3, logfile3, error=errlog, /get_lun ; logfile3 will be OUTPUT TABLE for elements by phase. openw, loglun4, logfile4, error=errlog, /get_lun IF (errlog NE 0) then PRINT, 'ERRROR opening logfile', !ERROR_STATE.MSG printf, loglun1, 'code: ', STRING(9B), codeName printf, loglun2, 'code: ', STRING(9B), codeName printf, loglun3, 'code: ', STRING(9B), codeName printf, loglun4, 'code: ', STRING(9B), codeName printf, loglun1, 'date: ', STRING(9B), SYSTIME() printf, loglun2, 'date: ', STRING(9B), SYSTIME() printf, loglun3, 'date: ', STRING(9B), SYSTIME() printf, loglun4, 'date: ', STRING(9B), SYSTIME() printf, loglun1, '================================' printf, loglun2, '================================' printf, loglun3, '================================' printf, loglun4, '================================' printf, loglun1, 'phases table: Pixels in each phase as identified by the code using current parameters.' printf, loglun1, FORMAT = '($%"%s \t %s \t %s \t %s \t %s \t")', 'name', 'total', 'mask', 'mask and holes', 'object' for x=1, (sizePhaseArray-1) do begin if cat[x].phsIncl then printf, loglun1, FORMAT = '($%"%s\t")', cat[x].label endfor for x=1, (sizePhaseArray-1) do begin if cat[x].phsIncl then printf, loglun1, FORMAT = '($%"%s\t")', cat[x].label endfor printf, loglun1, 'xx' printf, loglun1, FORMAT = '($%"%s \t %s \t %s \t %s \t %s \t")', '', 'pixels', 'pixels', 'pixels', 'pixels' for x=1, (sizePhaseArray-1) do begin if cat[x].phsIncl then printf, loglun1, FORMAT = '($%"%s\t")', 'pxlCount' endfor for x=1, (sizePhaseArray-1) do begin if cat[x].phsIncl then printf, loglun1, FORMAT = '($%"%s\t")', 'fraction' endfor printf, loglun1, 'xx' printf, loglun2, 'elsTotal table: min, max, stdev values for elements over all pixels in each object, and total intensity INCLUDING holes and unknowns.' printf, loglun2, FORMAT = '($%"%s \t %s \t %s \t %s \t %s \t")', 'name', 'total', 'mask', 'mask and holes', 'object' for z=1, (sizeElmArray-1) do begin ;Mask index is 0. printf, loglun2, FORMAT = '($%"%s \t %s \t %s \t %s \t %s \t")', elms[z].label, elms[z].label, elms[z].label, elms[z].label, elms[z].label endfor printf, loglun2, 'kk' printf, loglun2, FORMAT = '($%"%s \t %s \t %s \t %s \t %s \t")', '', 'pixels', 'pixels', 'pixels', 'pixels' for z=1, (sizeElmArray-1) do begin ;Mask index is 0. printf, loglun2, FORMAT = '($%"%s \t %s \t %s \t %s \t %s \t")', 'min', 'max', 'mean', 'stdev', 'total' endfor printf, loglun2, 'kk' ; Header for loglun3: elements by phase. ONLY through element S. One line per object measured. printf, loglun3, 'elsByPha table: Total intensity of each element in each identified phase over all unmasked pixels in each object.' printf, loglun3, FORMAT = '($%"\t%s\t")', 'phases>' for x=indxUnk, (sizePhaseArray-1) do begin ; ASSUME indxMsk = 0, indxUnk = 1. if cat[x].phsIncl then begin for z=1, (indxNi-1) do begin ; was to (numImages-1) which includes all elements. Mask index is 0. printf, loglun3, FORMAT = '($%"%s\t")', cat[x].label endfor endif endfor printf, loglun3, 'xx' ; This forces going to next line. ASSUMPTION: Order of phases is mask, unknowns, holes = 0, 1, 2 indices. printf, loglun3, FORMAT = '($%"n\t%s\t")', 'object' for x=indxUnk, (sizePhaseArray-1) do begin if cat[x].phsIncl then begin for z=1, (indxNi-1) do begin ; was to (numImages-1) which includes all elements. Mask index is 0. printf, loglun3, FORMAT = '($%"%s\t")', elms[z].label endfor endif endfor printf, loglun3, 'xx' ; This forces going to next line. ; READ CSV FILE ======================================================= ; rem code Input basename Maskname targetName map size X map size Y dwell time (ms) current (nA) csvWDSsetup ; INDICES of variables in CSV file (column number minus 1): iCodeColumn1 = 0 iBaseFilename = 1 iMaskFilename = 2 iTargFilename = 3 ixPxlperMap = 4 iyPxlperMap = 5 iDwellTime = 6 iBeamCurr = 7 iWDSsetup = 8 ;iWDSsetup = 2 ; Open file for formatted input sample list: openr, csvlun, csvInputFile, error=error, /get_lun IF (error NE 0) then PRINT, !ERROR_STATE.MSG lineIn = '' nVarsInCsv = 9 csvOffsets = intarr(nVarsInCsv) csvStrLens = intarr(nVarsInCsv) numSamples = 0 pos = 0 ; Scroll through open file. Store info for each item. ; NEW ORDER: basefile, targetfile, DwellTime, beamCurrentnA ; format: nVarsInCsv = # of rows ; First line has ONLY integer # of samples. Make arrays for them. readf, csvlun, numSamples ; Declare arrays based on number of samples: csvfilenum = strarr(numSamples) csvCodeColm1 = strarr(numSamples) csvBasefile = strarr(numSamples) csvMaskfile = strarr(numSamples) csvTargetfile = strarr(numSamples) csvxPxlPerMap = intarr(numSamples) csvyPxlPerMap = intarr(numSamples) csvmsDwellTime = intarr(numSamples) csvbeamCurrent = intarr(numSamples) csvWDSsetup = intarr(numSamples) csvNaPresent = intarr(numSamples) ; csvCAInum = strarr(numSamples) ; This is NOT used (yet) ; fileSpec = strarr(numSamples) ; date = make_array(numSamples, /STRING, VALUE = 'XXX-XX') ; Read in and parse the CSV file. for n=0, numSamples-1 do begin ; Input whole line (lineIn) from csv file with identity 'csvlun' readf, csvlun, lineIn print, 'n=',n, ' Line in>', lineIn, '< len = ', strlen(lineIn) csvOffsets = strsplit(lineIn, ',',LENGTH=csvStrLens) ;strsplit is an IDL function. ; Check size of results: print, ' sample =',n, 'Line in>', lineIn, '< len = ', strlen(lineIn) s = size(csvOffsets) print, ' size csvOffsets = ', s[1] s = size(csvStrLens) print, ' size csvStrLens = ', s[1] ; set values from imput line: basefile, targetfile, WDSsetup, PxlPerMap, DwellTime csvCodeColm1[n] = strmid(lineIn, csvOffsets[iCodeColumn1], csvStrLens[iCodeColumn1]) csvBasefile[n] = strmid(lineIn, csvOffsets[iBaseFilename], csvStrLens[iBaseFilename]) csvMaskfile[n] = strmid(lineIn, csvOffsets[iMaskFilename], csvStrLens[iMaskFilename]) csvTargetfile[n] = strmid(lineIn, csvOffsets[iTargFilename], csvStrLens[iTargFilename]) csvxPxlPerMap[n] = strmid(lineIn, csvOffsets[ixPxlperMap], csvStrLens[ixPxlperMap]) csvyPxlPerMap[n] = strmid(lineIn, csvOffsets[iyPxlperMap], csvStrLens[iyPxlperMap]) csvmsDwellTime[n] = strmid(lineIn, csvOffsets[iDwellTime], csvStrLens[iDwellTime]) csvbeamCurrent[n] = strmid(lineIn, csvOffsets[iBeamCurr], csvStrLens[iBeamCurr]) csvWDSsetup[n] = strmid(lineIn, csvOffsets[iWDSsetup], csvStrLens[iWDSsetup]) csvfilenum[n] = n csvNaPresent[n] = 0 ; OLD if ( strmid(lineIn, csvOffsets[iNaPresent], csvStrLens[iNaPresent]) EQ 1) then csvNaPresent[n] = TRUE else csvNaPresent[n] = FALSE ; THIS assumes that the filename contains only ONE underscore!! ; This is NOT used (workaround was found) ;csvOffsets = strsplit(csvBasefile[n], '_',LENGTH=csvStrLens) ;strsplit is an IDL function. ;csvCAInum[n] = strmid(lineIn, csvOffsets[1], csvStrLens[1],/REVERSE_OFFSET) endfor ; arrays for elements and phases allelmSum = DBLARR(sizePhaseArray, numSamples, sizeElmArray) ;allelmSum[i,j,k] where i = phase, j = object, k = sum allelmMean= fltarr(sizePhaseArray, numSamples, sizeElmArray) ;allelmMean[i,j,k] where i = phase, j = object, k = mean ;objectnames = strarr(numsamples) ;eleinphase = fltarr(numsamples, numImages, numImages) ; print to log file: printf, loglun4, 'sample-n:', numSamples printf, loglun4, 'Sample listing:' ; 1 2 3 4 5 6 7 8 printf, loglun4, FORMAT = '(%"\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t")', "n", "baseFile", "maskfile", "target", "Xpxls/map", "dwell(ms)", "curr(nA)", "Na?" for n=0, numSamples-1 do begin ; 1 2 3 4 5 6 7 8 printf, loglun4, FORMAT = '(%"CSV\t%i\t%s\t%s\t%s\t%i\t%i\t%i\t%s\t")', $ n, csvBasefile[n], csvMaskfile[n], csvTargetfile[n], csvxPxlPerMap[n], csvmsDwellTime[n], csvbeamCurrent[n], csvNaPresent[n] endfor flush, loglun4 ; The CSV input file is no longer necessary. Delete it. free_lun, csvlun ; end of input of data from .csv file =================================================== ; ======================================================================================= ; Values that should be set from CSV file: ; dwellTime = 60.0 ; milliseconds (ms) (r1 = 30ms, 20nA) ; TRUNCATE the number of elements 20-Jun-2016. ; fileSpecs = ['xx', 'Mg', 'Ca', 'Al', 'Si', 'Ti', 'Fe', 'S', 'Ni', 'Na'] ; 10 'elements' (images) ; numFileSpecs = 10 fileSpecs = ['xx', 'Mg', 'Ca', 'Al', 'Si', 'Ti', 'Fe'] ; 7 'elements' (images) numFileSpecs = 7 ; Fixed order for elements & maps: >> 'xx', 'Mg', 'Ca', 'Al', 'Si', 'Ti', 'Fe', (not present: 'S', Ni') fileSpecEnd = '_' numImages = 0 FOR i= 0, (sizeElmArray-1) DO elms[i].present = 0 FOR i= 0, indxFe DO BEGIN ; TRUNCATED! was indxNa until 20-Jun-2016. elms[i].present = 1 numImages = numImages + 1 ENDFOR if numImages NE 7 then print, 'ERROR! Wrong # of elements !! numImages=', numImages noMask = FALSE ; This value (inputMaskGS) is necessary if the mask was not set up correctly. ; inputMaskGS = 5 MapCrit = 2 ; for 32-bit 60ms/pxl 20nA data. ; Flags, etc boxScoreyes1no2 = 2 pxlIncrement1 = 1 ; To speed up tests, change from 1 to a larger number. Will default to 1. ; NOTE: If background is white, then auto levels will not brighten in photoshop!! ;inputBitDepth = 32 ;whiteMosaicBackYes1no2 = 2 ;autoLevelsYes1No2 = 2 ;bytscaleOn1Off2 = 2 ; add 18-Dec-2011 - turn OFF for 32-bit work. ;output32bitElmapsyes1no2 = 1 ; add 26-Jan-2012 ;rawImageStatsYes1no2 = 1 ; Print input data for batch processing: print,'n, basefile, maskfile, targetfile, xpxlperMap, xpxlperMap, DwellTime, beamCurrentnA, elSetup' for n=0, numSamples-1 do begin print, n, ',', csvBasefile[n], ',', csvMaskfile[n], ',', csvTargetfile[n], ',', csvmsDwellTime[n], ',', csvbeamCurrent[n], ',', csvxPxlPerMap[n], ',', csvWDSsetup[n] endfor curSampleN = 0 endelse ; END of segment for batch input 'inputFromListYes1no2' ; *********************************************************************************************************** printf, loglun4, FORMAT = '(%"\tnorm dwell time (ms) =\t%f\t")', criterionDwellTime printf, loglun4, FORMAT = '(%"\tnorm beam curr (nA) =\t%f\t")', criterionBeamCurrent printf, loglun4, FORMAT = '($%"\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t")', '', 'targetFile', 'norm-nA', 'norm-ms', 'crop-X', 'crop-Y', 'dwell(ms)', 'curr(nA)' printf, loglun4, 'xx' flush, loglun4 StartNextFileFromList: if (inputFromListYes1no2 EQ 1) then begin while (curSampleN NE numSamples) do begin ; Zero values for counts of phase pixels for p=0, (sizePhaseArray-1) do cat[p].count = 0 ; check the sample number ;if (curSampleN EQ numSamples) then goto, EndRun curSampleN = curSampleN + 1 print, ' Starting sample number ', curSampleN ; sourceBaseDir is already set ; sourceSubDir is already set ; targetBaseDir is already set basefile = csvBasefile[curSampleN-1] targetfile = csvTargetfile[curSampleN-1] if (flagMaskDifferYes1no2 EQ 1) then maskfile = csvMaskfile[curSampleN-1] ; add 19-July-2016 DSE dwellTime = csvmsDwellTime[curSampleN-1] beamCurrent = csvbeamCurrent[curSampleN-1] xPxlPerMap = INTARR(1) xPxlPerMap[0] = csvxPxlPerMap[curSampleN-1] fileNum = [ csvfilenum[curSampleN-1] ] ; WDSsetup = csvWDSsetup[curSampleN-1] ; Items set based on inputs from the CSV file: yPxlPerMap = xPxlPerMap mapSizeX = xPxlPerMap mapSizeY = yPxlPerMap title = csvTargetfile[curSampleN-1] fileSpecEnd = '_' ; NEW coding for setup 25-July-2016 (DSE) elSetup = csvWDSsetup[curSampleN-1] print, 'Code for elSetup = ', elSetup CASE elSetup OF 1: BEGIN END 2: BEGIN ; For ECP standard WDS + EDS element map set fileSpecs = ['xx', 'Mg', 'Ca', 'Al', 'Si', 'Ti', 'Fe', 'S', 'Ni'] ; , 'Na'] ; Fixed order for elements & maps: >> 'xx', 'Mg', 'Ca', 'Al', 'Si', 'Ti', 'Fe', (not present: 'S', Ni') fileSpecEnd = '_' numImages = 0 FOR i= 0, (sizeElmArray-1) DO elms[i].present = 0 elms[indxXx].present = 1 elms[indxMg].present = 1 elms[indxCa].present = 1 elms[indxAl].present = 1 elms[indxSi].present = 1 elms[indxTi].present = 1 elms[indxFe].present = 1 elms[indxS].present = 1 elms[indxNi].present = 1 numImages = 9 6: BEGIN ; For ECP standard WDS + EDS element map set fileSpecs = ['xx', 'Mg', 'Ca', 'Al', 'Si', 'Ti', 'Fe', 'S', 'Ni', 'Na'] fileSpecEnd = '_' numImages = 0 FOR i= 0, (sizeElmArray-1) DO elms[i].present = 0 elms[indxXx].present = 1 elms[indxMg].present = 1 elms[indxCa].present = 1 elms[indxAl].present = 1 elms[indxSi].present = 1 elms[indxTi].present = 1 elms[indxFe].present = 1 elms[indxS].present = 1 elms[indxNi].present = 1 elms[indxNa].present = 1 numImages = 10 END ELSE: print, 'elSetup has illegal value.' ENDCASE print, '---------------' MapCrit = 2 CASE MapCrit OF 0: BEGIN print, ' SKIP map criteria ....' for p=0, (sizePhaseArray-1) do cat[p].phsIncl = FALSE END ; REMOVED 27aug2012: all MapCrit EQ 1 mode determinations from 22-Aug-2012 codes. 2: BEGIN ; This is the NEW scheme. ECP modified October 2016. DSE starts this 2-Aug-2012. For 32-bit 60ms/pxl 20nA data. print, 'FOR CAIs!' mapCriterion = 'CAIs-scheme1' criterionDwellTime = 60.0 criterionBeamCurrent = 20.0 criterionKeV = 15 for p=1, (sizePhaseArray-1) do cat[p].phsIncl = FALSE ; ECP phase list: Spn/TO2/Prv/Ttn/SiO2/Wol/Cor/Hib/Olv/Opx/Cpx/Mel/Plg/Ilm/Ab cat[indxMsk ].phsIncl = TRUE cat[indxUnk ].phsIncl = TRUE cat[indxHole ].phsIncl = TRUE ; Set holes = mask color-1 (change from mask color 12-July-2016 DSE). cat[indxCc ].phsIncl = FALSE cat[indxMet ].phsIncl = TRUE cat[indxFeS ].phsIncl = TRUE cat[indxCaS ].phsIncl = FALSE cat[indxMnS ].phsIncl = FALSE cat[indxMgS ].phsIncl = FALSE cat[indxSpn ].phsIncl = TRUE cat[indxTO2 ].phsIncl = TRUE ;TiO2 cat[indxPrv ].phsIncl = TRUE ;CaTiO3 cat[indxTtn ].phsIncl = FALSE ;sphene-CaTiSiO5 = titanite cat[indxSiO2 ].phsIncl = TRUE cat[indxWol ].phsIncl = FALSE cat[indxCor ].phsIncl = TRUE cat[indxHib ].phsIncl = TRUE cat[indxCA2 ].phsIncl = FALSE cat[indxCA1 ].phsIncl = FALSE cat[indxOlv ].phsIncl = TRUE cat[indxGls ].phsIncl = TRUE cat[indxOpx ].phsIncl = TRUE cat[indxCpx ].phsIncl = TRUE cat[indxMel ].phsIncl = TRUE cat[indxPlg ].phsIncl = TRUE cat[indxAlt ].phsIncl = FALSE cat[indxRhn ].phsIncl = FALSE cat[indxIlm ].phsIncl = TRUE cat[indxAb ].phsIncl = TRUE cat[indxHrc ].phsIncl = FALSE cat[indxCaTs ].phsIncl = FALSE cat[indxMsc ].phsIncl = FALSE ; Criteria October 2016 (modified from 2012) ;Holes and cracks holeMinAll = 400 ;400.0 works ok but misses some ; need this! ;holeMinAll = 100.0 ; for standard grid ; sulfide S_AlSi (S - Al - Si) (old SminSulfide) sulMinS = 100.0 ; was 240 from standards ; metal F_CMAS (Fe - sum of Ca+Mg+Al+Si > this number) (old FmCMASminMetal) metMinFe = 100.0 ; quartz minimum in Si - 2Mg-Fe-Al-2Ca-Ti (Sim2MgFeAl2CaTi map) SimSiO2 = 700.0 ; TiO2 minimum in Ti - 2Mg-2Fe-2Al-2Ca-2Si (Tim2MgFeAlCaSi map) TinTiO2 = 1400.0 ; corundum minimum in Al - 2Si-2Mg-Fe-2Ca (Alm2Si2MgFe2Ca map) Alm2S2MF2C = 850 ;1000 works well for samples ; wollastonite CaSiO3 minCSwoll = 630 ; in Ca + Si - 3*(Al+Mg+Fe) (silica is higher) ; perovskite CaTiO3 ; wasCa3Ti3m5Al2Si2FePRV = 175 Was 200 until 29-Aug-2012. 800 works better for stds when sphene + ilmenite present. prvMinCaTi = 100.0 ; 200; 1100 for standard grid 32-bit value for 60ms 20nA in Ca3Tim5Al2Si2FePRV map (8bit = 42) ; Spinel MgAl2O4 in (Mg+Al) - (5*Ca+2*Si) map (from 9/12/2011) ;MgAlm2Si5CaminSpn = 150.0 ; 650 for standard grid 32-bit value for 60ms 20nA in MgAlm2Si5Ca map (8bit = 6) ;CM2SoverAminSpn = 0.02 ; 0.09 for standard grid Ca+Mg-2Si over Al maxCaspn = 19.9 AMm2SminSpn = 250 ; 250 works well Al+Mg-2Si AMm2SmaxSpn = 1200 ;950 is good for samples ;730 ;AMm2SoverCminSpn = 76.0 ;Al+Mg-2Si over Ca ;AMm2SoverCminSpn = 76.0 ;standard grid ; Sphene CaTiSiO5 Ca+Ti+Si - 5*(Al+Mg+Fe) minCTSsphene = 890.0 minTisphene = 600.0 ; Ilmenite FeTiO3 Fe+Ti - 2*(Si+Ca+Al) minFTilmenite = 850;0 minTiIlm = 800.0 ; carbonate criteria in C_ASMF (Ca - Al - Si - Mg - Fe) ccMinCa = 95 ; Ca-Al-Ti only phases; Hib, CA1, CA2 minCATmMSF = 500.0 ; This is not going to work for CA1 and CA2... we have no standards for them! CAmMminHib = 600 ;AloCaminHib = 4.0 ; minimum Al/Ca (range = 15 - 30) ;AloCamaxHib = 35.0 AloCaminCA2 = 2.0 AloCaminCA1 = 0.1 ; hercynite (Al+Fe-2Si)/Ca AFm2SoverCminHrc = 0.6 ; olvine and opx MgFeSim3Ca2AlTiminOLOPX = 300.0 ; 340 for standard grid ;MgFeSim3Ca2AlTiminOLOPX = 200.0 ; standard grid MgFevsSiPX = 0.4 ; in ratio map: MgFeoverSi.tif MgFevsSiOLmin = 1.5 ;1.5 good for samples and 1.6 for stdgrids, was 1.4 MgFevsSiOLmax = 5 ;4 worked well in samples ;MgFevsSiOLmin = 0.5 ;standard grid ; feldspar minSiplg = 0.5 ;minSiplg = 500 ;standard grid minplag = 5.0 ;2.0 for standard grid minimum SAover2CM for plagioclase, this separates plag from mel and cpx ;minplag = 2.0 ;standard grid ; DSE Part for plagioclase (anorthite) and other feldspar. minSCAplg = 1200 ;1500.0 works well for samples ;was 1700.0 minSoAplag = 0.1 maxSoAplag = 0.2 maxSoAalb = 0.45 minNaAlbite = 50 ; was 130.0 ; Cpx + Mel + An + mesostatsis minCaCpxMelAn = 20.0 ;80.0 when including Fe minCameso = 10 ; 80 for standard grid minimum Ca in these minerals. In Ca - S map if S present. ;CaMinCoMASFmel = 0.80 ; for samples min allowed Ca/(Mg+Al+Si+Fe) in mel (max 0.6) ;CaMinCoMASFmel = 0.45 ;for standard grid ;CoMSminmel = 1.0 ; for standard grid ;CoMSmincpx = 0.15 ;SoCminmel = SoCmaxmel = 1.59 ;was 1.7 SoCminmel = 0.1 ;SAMCmaxmel = 120 SoCmincpx = 1.6 ;was 1.8 SoCmaxcpx = 6.0 ;was 4.0, 5.0 worked well ;SAMCmaxcpx = 265 ;SAMCmincpx = 100 SACmaxcpx = 400 ;400 ;MSCmincpx = 430 ;to distinguish between cpx and mesostasis in chondrules SoCminmeso = 4.0 ;was 5.01, 6.0 SoCmaxmeso = 60 ;CoMSmincpx = 0.5 ; for standard grid ;CaMinCoMASFcpx = 0.09 ;0.13;30 ; 0.19 for standard grid min allowed Ca/(Mg+Al+Si+Fe) in diopside augite.is lower ;CaMinCoMASFcpx = 0.5 ; 0.19 with A and F for standard grid ;CaMinCoMASFplg = 0.009 ; 0.15 for standard grid max allowed Ca/(Mg+Al+Si+Fe) in plag. (max 0.19) ;CaMinCoMASFplg = 0.15 ; for standard grid ;CoMSminplg = 0.0009 CoMSminplg = 0.015 ; for standard grid maxFePlg = 9.0 ; max allowed Fe in plag ; Ca-Al rich pyroxene ACoverSmaxCaTs = 2.6 ACoverSminCaTs = 2.0 ; Make ratio maps. Add a constant to each denominator to avoid division by zero errors. ratConstDenom = 0.01 END ELSE: print, 'MapCrit has illegal value!' ENDCASE ; Check for obvious errors in criteria: ;if (beamCurrent EQ criterionBeamCurrent) then normBeamCurrent = FALSE else normBeamCurrent = TRUE ;if (dwellTime EQ criteriondwellTime) then normDwellTime = FALSE else normDwellTime = TRUE print, ' ------' ; Declare an array of unsigned long integers to store nElements by 3 (x,y,val) values. maxXY = ULONARR(numImages, 3) print, 'grayscale totals for elements: ', sourceBaseDir print, 'file elem NpxlTot normVal nWht nBlk nGray totVal' ; Load all maps **************************************************************************************************** histoNbins = 100 histoMax = 1000 histo = FLTARR(numImages, histoNbins) ; fileSuffix = '' for n=0, (numImages-1) do begin IF (elms[n].present) THEN BEGIN if ( (flagMaskDifferYes1no2 EQ 1) AND (n EQ 0) ) then begin ; Add 19-July-2016 DSE ; --- Construct the name of the MASK file. nameOnly = maskfile + '.tif' ; csvMaskfile[n] name = maskFileDir + nameOnly endif else begin ; if a mask is present, but in a separate directory, enter it here. (20jun2016) if (noMask EQ FALSE) then begin if (n EQ 0) then begin ; --- Construct the name of the MASK file. nameOnly = basefile + fileSpecEnd + fileSpecs[n] + '.tif'; fileSuffix + name = maskFileDir + nameOnly endif else begin ; --- Construct the name of the ELEMENT file. nameOnly = basefile + fileSpecEnd + fileSpecs[n] + '.tif'; fileSuffix + name = sourceELMmapSubDir + nameOnly endelse endif ; if NO mask is present, skip n=0, and build mask layer later! (26feb2012) if ( (noMask EQ TRUE) AND (n EQ 0) ) then n = 1 endelse ; --- Input the TIFF data. 32-bit EMP files are FLOAT type (type 4) floating point #s. print, 'Reading file ', name, ' Element = ', elms[n].label, ' present = ', elms[n].present slice = READ_TIFF(name) ; --- Obtain size of the map. maxXY[n, 2] = MAX(slice) ; get the size (number of pixels) s = SIZE(slice) print, ' size(0)=dims=',s[0], ' dim(1)=',s[1],' dim(2)=', s[2], ' (byte=0, int=2, longint=3)', ' type(3)=',s[3], ' tot size(4)=',s[4] ;print, ' type(3)=',s[3], ' tot size(4)=',s[4] IF ( (n EQ 0) OR ((noMask EQ TRUE) AND (n EQ 1)) ) THEN BEGIN maxElmapX = s[1] maxElmapY = s[2] ; declare arrays of the correct size: (25-Feb-2012: Change xmaps from ULONARR to FLTARR) Changed again 8-Aug-2012. if (inBitDepth NE 8) AND (inBitDepth NE 32) then print, 'ERROR! inBitDepth NOT 8 or 32 !!!!' & WAIT, 10.0 if (inBitDepth EQ 8) then xmaps = UINTARR(numImages, maxElmapX, maxElmapY) if (inBitDepth EQ 32) then xmaps = FLTARR(numImages, maxElmapX, maxElmapY) ENDIF ; --- Place map in multi-element array. ; To preserve 32-bit input depth, use 'ULONG' to make these 32-bit. Otherwise, use 'UINT'. if (inBitDepth EQ 8) then xmaps[n,*,*] = UINT(slice[*,*]) if (inBitDepth EQ 32) then xmaps[n,*,*] = FLOAT(slice[*,*]) ; 32-bit unsigned integer array type = ULONG. 16-bit uns. int = UINT. Raw maps are float type. ; --- (optional, debug) send this map to screen. ; TV, xmaps[n,*,*], /ORDER ; print, 'name...', name goto, bypassHisto ; View a histogram for this map: print, 'Histogram of file ', fileSpecs[n], ' MIN = ', MIN(xmaps[n,*,*]), ' MAX = ', MAX(xmaps[n,*,*]) histo(n, *) = histogram(xmaps[n, *, *], NBINS=histoNbins, MAX=histoMax) sHisto = SIZE(histo) print, 'HISTO: dims=',sHisto[0], ' dim(1)=',sHisto[1],' dim(2)=', sHisto[2], ' type =',sHisto[3], ' tot size =',sHisto[4] print, ' histo(n, 0) = ', histo(n, 0) PLOT, histo(n, *) ;, MAX_VALUE=1000.0 print, '---------------' bypassHisto: goto, bypassDebug1 ; debugging output 8aug2012: CTSm5AMF if (n GE indxFe) then begin if (inBitDepth EQ 8) then r = INTARR(1, maxElmapX, maxElmapY) else r = FLTARR(1, maxElmapX, maxElmapY) r = xmaps[indxCa,*,*] + xmaps[indxTi,*,*] + xmaps[indxSi,*,*] + (-5.0) * xmaps[indxAl,*,*] + (-5.0) * xmaps[indxMg,*,*] + (-5.0) * xmaps[indxFe,*,*] r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = targetBaseDir + targetfile + 'TEST_1_CTSm5AMF.tif' if (inBitDepth EQ 8) then begin r[0,*,*] = BYTSCL(r[0,*,*]) if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse print, ' ----------' endif bypassDebug1: ; goto, bypassNorms ; Normalize to dwell time 120ms. But NOT for the mask! if ( normDwellTime AND (n NE indxXx) ) then begin if (dwelltime EQ 0.0) then print, 'ERROR! dwelltime is zero! STOP STOP ERROR ERROR ERROR' if (criterionDwellTime EQ 0.0) then print, 'ERROR! criterionDwellTime is zero! STOP STOP ERROR ERROR ERROR' factor = FLOAT(criterionDwellTime/FLOAT(dwellTime)) print, ' element map ', elms[n].label, ' normalized for criterion dwell time = ', criterionDwellTime, ' relative to actual dwell time = ', dwellTime, $ ' Beam current normalization factor = ', factor, ' (criterionDwellTime/dwellTime).' xmaps[n,*,*] = factor * xmaps[n,*,*] print, '-----' endif ; Normalize to beam current. But NOT for the mask! if ( normBeamCurrent AND (n NE indxXx) ) then begin if (beamCurrent EQ 0.0) then print, 'ERROR! beamCurrent is zero! STOP STOP ERROR ERROR ERROR' if (criterionBeamCurrent EQ 0.0) then print, 'ERROR! criterionBeamCurrent is zero! STOP STOP ERROR ERROR ERROR' factor = (criterionBeamCurrent/beamCurrent) print, ' element map ', elms[n].label, ' normalized for criterion beam current = ', criterionBeamCurrent, ' relative to actual current = ', beamCurrent, $ ' Beam current normalization factor = ', factor, ' (criterionBeamCurrent/beamCurrent).' xmaps[n,*,*] = xmaps[n,*,*] * (criterionBeamCurrent/beamCurrent) print, '-----' endif ENDIF ; End of conditional 'if elms[n].present' ; bypassNorms: endfor ; End of loop on images (maps) print, ' output array size is dim(1)=x=', maxElmapX,' dim(2)=y=', maxElmapY goto, bypassDebug2 ; debugging output 8aug2012: CTSm5AMF if (n GE indxFe) then begin if (inBitDepth EQ 8) then r = INTARR(1, maxElmapX, maxElmapY) else r = FLTARR(1, maxElmapX, maxElmapY) r = xmaps[indxCa,*,*] + xmaps[indxTi,*,*] + xmaps[indxSi,*,*] + (-5.0) * xmaps[indxAl,*,*] + (-5.0) * xmaps[indxMg,*,*] + (-5.0) * xmaps[indxFe,*,*] ; r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = targetBaseDir + targetfile + 'TEST_2_CTSm5AMF.tif' if (inBitDepth EQ 8) then begin r[0,*,*] = BYTSCL(r[0,*,*]) if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse print, ' ----------' endif bypassDebug2: ; If NO mask, build that layer now. Make all the mask 'unknown' if (noMask EQ TRUE) then begin xmaps[indxXx,*,*] = gsUnk print, 'Histogram of file ', fileSpecs[indxMsk], ' MIN = ', MIN(xmaps[indxMsk,*,*]), ' MAX = ', MAX(xmaps[indxMsk,*,*]) histo(indxXx, *) = histogram(xmaps[indxMsk, *, *], NBINS=histoNbins, MAX=histoMax) endif goto, bypassHistoLog ; Output histograms to log file: printf, loglun, 'Histogram Nbins = ', histoNbins, ' Max = ', histoMax for n=0, numImages-1 do printf, loglun, FORMAT = '($%"%s \t")', fileSpecs[n] printf, loglun, 'xx' for i=0, histoNbins-1 do begin for n=0, numImages-1 do begin printf, loglun, FORMAT = '($%"%f \t")', histo[n, i] endfor printf, loglun, 'xx' endfor bypassHistoLog: print, '-----------' ; ************************************************************ ; auto_levels: This is legacy code not used here. ;calls the procedure stored at path: E;\tomo1\IDLprogs\auto_Level_v1.pro ;PRO auto_Level_v1, elMap, minCutoff, maxCutoff < arguments ; NOTE: autoLevels does NOT work if whiteMosaicBack is ON!! ;maxCutoff = 0.995 ;minCutoff = 0.995 ;auto_Level_v2, maskMap, minCutoff, maxCutoff ;print, ' returned from autoLevelproc' ; ********************************************************************************************************* ; Count the number of phases included. nPhsIncl = 0 for x=0, (sizePhaseArray-1) do if (cat[x].phsIncl EQ 1) then nPhsIncl = nPhsIncl + 1 if (boxScoreyes1no2 EQ 1) then begin ; ALTERNATE way to get numbers: ; step 1: Count in the 'box' for selected elements and phases: ; boxCnt will have min, max, mean, variance, stdev, as 0-5 items for each element. boxCnt = FLTARR(sizePhaseArray, numImages, 5) for p=0, (sizePhaseArray-1) do begin if (box[p, 0] GT 1) then begin for i=0, (numImages-1) do begin if (elms[i].present EQ 1) then begin ; loop through the x-y values in the box, and add up pixels. print, 'Counting el ', elms[i].label, ' in box: x = ', box[p, 0], ' to ', box[p, 1], ' y = ', box[p, 2], ' to ', box[p, 3] boxCnt[p, i, 0] = MIN(xmaps[i, box[p, 0]:box[p, 1], box[p, 2]:box[p, 3]]) boxCnt[p, i, 1] = MAX(xmaps[i, box[p, 0]:box[p, 1], box[p, 2]:box[p, 3]]) result = MOMENT(xmaps[i, box[p, 0]:box[p, 1], box[p, 2]:box[p, 3]]) boxCnt[p, i, 4] = STDDEV(xmaps[i, box[p, 0]:box[p, 1], box[p, 2]:box[p, 3]]) boxCnt[p, i, 2] = result[0] boxCnt[p, i, 3] = result[1] endif print, ' box counted for element ', elms[i].label, ' MIN-MAX= ', boxCnt[p, 0], ' - ', boxCnt[p, 1] endfor ; end of loop on element maps ; set limits from this data: REMOVED! endif endfor print, 'elm', STRING(9B), 'MIN', STRING(9B), 'MIN', STRING(9B), 'AVG', STRING(9B), 'StdDev', STRING(9B), 'Variance' for p=0, (sizePhaseArray-1) do begin if (box[p, 0] GT 1) then begin print, cat[p].label, ' ======================' for i=0, (numImages-1) do begin if (elms[i].present EQ 1) then begin if (box[p, 0] GT 0) then begin print, elms[i].label, STRING(9B), boxCnt[p, i, 0], STRING(9B), boxCnt[p, i, 1], STRING(9B), $ boxCnt[p, i, 2], STRING(9B), boxCnt[p, i, 4], STRING(9B), boxCnt[p, i, 3] endif endif endfor endif endfor print, '----------------------------' endif ; end of conditional 'boxScoreyes1no2' ; ================================================================================================================== ; PRELIMINARY stuff for ALL 'MapCrit' choices: print, 'Start loop on all pixels for preliminary stuff.' ; NOT put in separate code: prelim_mask_set_v1 for x=0, (maxElmapX-1) do begin for y=0, (maxElmapY-1) do begin ; Set basic unknown and mask values. This ASSUMES that Mask is grayscale ZERO. if (xmaps[indxXx,x,y] EQ 0) then begin ; If this is a masked off pixel (black) in original mask, then don't go through the whole phase array. xmaps[indxXx,x,y] = cat[indxMsk].gsVal1 endif else begin ; If this is not a masked off pixel, make it an unknown in the result. xmaps[indxXx,x,y] = gsUnk endelse endfor ; end of loop on y (all pixels) endfor ; end of loop on x (all pixels) print, 'Finished preliminary stuff. count[indxMsk] = ', cat[indxMsk].count ; DEBUG ;map = xmaps[indxXx, *, *] ;TV, map ;print, ' Done with setting masked pixels ... ' ; ================================================================================================================== ; Find the dimensions of the masked area ************************************************************************** if (cropMapsyes1no2 EQ 1) then begin print, ' xmaps[indxXx] MASK dims = ', maxElmapX, ' by ', maxElmapY maskMap = FLTARR(maxElmapX, maxElmapY) maskMap[*,*] = xmaps[indxXx,*,*] xyLim = INTARR(2, 2) ;xyLim[2,2]: cols: x, y rows: min, max. print, 'Calling external locate_object_xylims_v1 routine to get xyLim values... ' ; CALL SUBPROGRAM locate_object_xylims_v1, maskMap, cat[indxMsk].gsVal1, xyLim maskMap = 0 print, 'Found x-y limits = xMin= ', xyLim[0,0], ' xMax= ', xyLim[0,1], ' yMin= ', xyLim[1,0], ' yMax= ', xyLim[1,1] print, '---------------' endif ; CROP MAPS ======================================================================================================== ; Crop to the dimensions of the masked area + buffer pixels each side *************************************************************************************** if (cropMapsyes1no2 EQ 1) then begin ; was this call: crop_to_mask_v1, xmaps, xyLim, xyBuffer xyBuffer = cropBufferMapPxls ; usually set to 5 print, 'x-y limits = xMin= ', xyLim[0,0], ' xMax= ', xyLim[0,1], ' yMin= ', xyLim[1,0], ' yMax= ', xyLim[1,1] xyLim[0,0] = xyLim[0,0] - xyBuffer if (xyLim[0,0] LT 0) then xyLim[0,0] = 0 xyLim[0,1] = xyLim[0,1] + xyBuffer if (xyLim[0,1] GT (maxElmapX-1)) then xyLim[0,1] = maxElmapX-1 xyLim[1,0] = xyLim[1,0] - xyBuffer if (xyLim[1,0] LT 0) then xyLim[1,0] = 0 xyLim[1,1] = xyLim[1,1] + xyBuffer if (xyLim[1,1] GT (maxElmapY-1)) then xyLim[1,1] = maxElmapY-1 print, 'NEW x-y limits = xMin= ', xyLim[0,0], ' xMax= ', xyLim[0,1], ' yMin= ', xyLim[1,0], ' yMax= ', xyLim[1,1] szNewX = xyLim[0,1]-xyLim[0,0]+1 szNewY = xyLim[1,1]-xyLim[1,0]+1 print, 'NEW x-y size X: ', szNewX, ' Y: ', szNewY newMaps = FLTARR(numImages, szNewX, szNewY) for m=0, (sizeElmArray-1) do begin if (elms[m].present) then begin newMaps[m, 0:(szNewX-1), 0:(szNewY-1)] = xMaps[m, xyLim[0,0]:xyLim[0,1], xyLim[1,0]:xyLim[1,1]] ; TV, newMaps[m,*,*], /ORDER print, 'Done resetting map image for m= ', m, ' label = ', elms[m].label endif endfor print, 'Setting all map images equal to new maps. Also setting new size constraints. ' xMaps = newMaps maxCropElX = szNewX ; WAS 'maxElX = szNewX' prior to 30jun2016 maxCropElY = szNewY ; WAS 'maxElY = szNewY' prior to 30jun2016 newMaps = 0.0 ; DEBUG: ;s = SIZE(xMaps) ;print, 'after cropping xMaps: dimensions = ', s[0], ' size x = ', s[1], ' size y = ', s[2] ;print, ' typeCode = ', s[3], ' #elements = ', s[4] endif ; end of part to crop to the dimensions of the masked area ; ================================================================================================================== ; Eliminate Bright Pixels *************************************************************************************** ; Pixel-by-Pixel interrogation of maps! if (delBrightPxlyes1no2 EQ 1) then begin elMap = FLTARR(maxCropElX, maxCropElY) ; Assumes that indxXx = 0, that the mask is index m = 0 !!!!! ; While the number of total possible elements is sizeElmArray, the actual number present is numImages. for m=1, (sizeElmArray-1) do begin if (elms[m].present) then begin elMap[*,*] = xmaps[m,*,*] ; CALL SUBPROGRAM del_bright_pxls_v3, elMap, brightDiffCrit, brightThreshold, hotPxlFound xmaps[m,*,*] = elMap[*,*] ; TV, xmaps[m,*,*], /ORDER endif print, 'Finished eliminating bright pixels for m = ', m, ' label = ', elms[m].label if (hotPxlFound) then print, ' Hot (over-bright) pixel found in element map ', elms[m].label, ' !!!!' endfor elMap = 0 print, 'Finished eliminating bright pixels.' endif ; end of conditional on 'delBrightPxlyes1no2'. ; end of Eliminate Bright Pixels ; ================================================================================================================== ; Find Forsterite for calibration ================================================================================== ; greate new map... ; Part for OLIVINE and OPX ************************************************** IF (findForstYes1no2) THEN BEGIN ; Find Mg+Si rich pixels in if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxMg,*,*] + xmaps[indxSi,*,*] + (-3.0) * xmaps[indxCa,*,*] + (-3.0) * xmaps[indxAl,*,*] if (elms[indxFe].present) then r = r - (3.0) * xmaps[indxFe,*,*] if (elms[indxTi].present) then r = r - (3.0) * xmaps[indxTi,*,*] if (elms[indxS].present) then r = r - (3.0) * xmaps[indxS,*,*] ; set values of MgSim3CaFeAlTi less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'fo-MgSim3CaFeAlTi.tif' if (inBitDepth EQ 8) then begin r[0,*,*] = BYTSCL(r[0,*,*]) ; Scale r to 0-255: ;TV, r[0,*,*], /ORDER ;print, '...'s if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse r2 = FLTARR(1, maxCropElX, maxCropElY) ; Make ratio map. Add a constant to each denominator to avoid division by zero errors. (ratConstDenom = 0.01 set above) r2 = ( xmaps[indxMg,*,*] + xmaps[indxFe,*,*] ) / ( xmaps[indxSi,*,*] + ratConstDenom ) r2[0,*,*] = (r[0,*,*] GE MgFeSim3Ca2AlTiminOLOPX) * r2[0,*,*] outFile = CritMapDir + targetfile + 'MgFeoverSi.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r2[0,*,*], /FLOAT for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin ; Find Mg-Fe-Si rich pixels and then use ratio (Mg+Fe)/Si to tell OL and OPX apart. if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] GE MgFeSim3Ca2AlTiminOLOPX) then begin if (r2[0,x,y] GE MgFevsSiOLmin) AND (cat[indxOlv].phsIncl EQ TRUE) then begin xmaps[indxXx,x,y] = cat[indxOlv].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxOlv].rgbX[z] cat[indxOlv].count = cat[indxOlv].count + 1 endif else begin if (cat[indxOpx].phsIncl EQ TRUE) then begin xmaps[indxXx,x,y] = cat[indxOpx].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxOpx].rgbX[z] cat[indxOpx].count = cat[indxOpx].count + 1 endif endelse endif endif endfor ; end loop on y. endfor ; end loop on x r2 = 0 ; destroy r2 array. ENDIF ; ================================================================================================================== ; APPLY MASK ====================================================================================================== ; =========== option: Print 32-bit TIF of each element map, and the mask (mask will be black in these outputs), ; Assumes that indxXx = 0, that the mask is index m = 0 !!!!! ; While the number of total possible elements is sizeElmArray, the actual number present is numImages. for m=0, (sizeElmArray-1) do begin if (elms[m].present) then begin outFileSpec = maskedFileDir + title + '-' + elms[m].label + '-new.tif' ; Impose the mask ONLY on the element maps (not on the mask itself!) ; Here, the conditional is used to control the new map. if (m NE 0) then xmaps[m,*,*] = xmaps[m,*,*] * (xmaps[indxXx,*,*] NE cat[indxMsk].gsVal1) ; print, 'targetBaseDir = ', targetBaseDir, ' title=', title, ' elm=', elms[m].label if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFileSpec, xmaps[m,*,*], /FLOAT ;DEBUG ;TV, xmaps[m,*,*], /ORDER ;print, ' ... TV printed map with mask cover... for ', elms[m].label endif endfor outFileSpec = maskedFileDir + title + '-' + 'mask-8bt.tif' elmap = BYTARR(maxCropElX, maxCropElY) elmap[*,*] = xmaps[indxXx,*,*] if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFileSpec, elmap, 1 elmap = 0 ; done with this variable print, 'Finished printing cropped element maps with eliminated bright pixels and with black mask.' ; end of print segment for new element maps. ; Count the masked area: Pixel - By - Pixel =============================================================== ; Use elemAvgVal to collect info for element maps, cropped, and with black masks. ; NOTE: This gives WRONG values for MIN, MEAN and STDEV since the masked area is also counted. for n=0, (numImages-1) do begin elemAvgVal[n, 0] = 1000.0 ;MIN(xmaps[n,*,*]) dummy value for below elemAvgVal[n, 1] = MAX(xmaps[n,*,*]) elemAvgVal[n, 2] = 0.0 ;MEAN(xmaps[n,*,*]) elemAvgVal[n, 3] = 0.0 ;STDDEV(xmaps[n,*,*]) elemAvgVal[n, 4] = TOTAL(xmaps[n,*,*]) ; This will be ok since the mask is zero valued. elemAvgVal[n, 5] = DOUBLE(0.0) ; add 12-July-2016 DSE to count up totals later. endfor ; Count the number of pixels in the masked area: if (cropMapsyes1no2 EQ 1) then begin countTot = LONG(maxCropElX)*LONG(maxCropElY) ; add 3july2016 DSE. countTot must be a LONG integer! countObj = 0L for m=0, (maxCropElX-1) do begin for n=0, (maxCropElY-1) do begin if (xMaps[indxXx, m, n] EQ cat[indxUnk].gsVal1) then begin countObj = countObj + 1 for i=1, (numImages-1) do begin if (xMaps[i, m, n] LE elemAvgVal[i, 0]) then elemAvgVal[i, 0] = xMaps[i, m, n] endfor endif endfor endfor cat[indxMsk].count = countTot - countObj for i=1, (numImages-1) do elemAvgVal[i, 2] = elemAvgVal[i, 4]/countObj endif ; MODE CALCULATION: Pixel - By - Pixel ========================================================================== ; REMOVED: all MapCrit EQ 1 mode determinations from 22-Aug-2012 codes. ; *********************************************************************RRRR*************************** ; *********************************RRR*****************************S********************************** ; **********************RRR************************************************************************** ; ************************************************************************************************ ; ************************************************************************************************ ; ************************************************************************************************ ; ************************************************************************************************ ; **********************************************BYPASS MODE MAPPING************************************** if (bypassModeMappingYes1no2 EQ 1) then goto, bypassModeMapping print, ' Entering pixel-by-pixel counting part for map set #', curSampleN ; The arrays for elemSums and elemMeans have size 10. for p=0, (sizePhaseArray-1) do begin cat[p].count = 0 ; add 3july2016 DSE; make DBL 12july2016 for z=0, 9 do cat[p].elemSums[z] = DOUBLE(0.0) endfor ; To run tests, do in increments of 5 pixels: for x=0, (maxCropElX-1), 5 do begins if (pxlIncrement1 EQ 0) then pxlIncrement1 = 1 ; DEBUG: for p=1, (sizePhaseArray-1) do IF (cat[p].phsIncl EQ TRUE) THEN print, 'phase ', cat[p].label, ' PRESENT in counting.' print, '------' ; Set up modeRGB for final colored mode map: modeRGB = BYTARR(3, maxCropElX, maxCropElY) cat[indxMsk].count = 0 ; added 3-July-2016 DSE for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if (xmaps[indxXx,x,y] EQ cat[indxMsk].gsVal1) then begin for z=0, 2 do modeRGB[z, x, y] = cat[indxMsk].rgbX[z] cat[indxMsk].count = cat[indxMsk].count + 1 endif else begin ; If this is not a masked off pixel, make it an unknown in the result. for z=0, 2 do modeRGB[z, x, y] = cat[indxUnk].rgbX[z] endelse endfor ; end of loop on y (all pixels) endfor ; end of loop on x (all pixels) ; TV, modeRGB print, '...' print, 'Finished setup of modeRGB array map for final mode image. Masked pxls = ', cat[indxMsk].count ; REMOVED: all MapCrit EQ 1 map setups from 22-Aug-2012 codes. IF MapCrit EQ 2 THEN BEGIN print, ' ENTERING mode criteria #2.' ; Part for HOLES in the material, or cracks. Set holes = their own colors (changed 12-July-2016 DSE). ************************************************** IF (cat[indxHole].phsIncl EQ TRUE) THEN BEGIN ; Find pixels with very low counts if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxMg,*,*] + xmaps[indxFe,*,*] + xmaps[indxAl,*,*] + xmaps[indxSi,*,*] + xmaps[indxCa,*,*] + xmaps[indxTi,*,*] if (elms[indxS].present) then r = r + xmaps[indxS,*,*] if (elms[indxNi].present) then r = r + xmaps[indxNi,*,*] ; set values of sumAllEls less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'sumAllEls.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; Holes set to value of Mask. for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] LE holeMinAll) then begin xmaps[indxXx,x,y] = cat[indxHole].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxHole].rgbX[z] cat[indxHole].count = cat[indxHole].count + 1 endif endif endfor endfor ENDIF ; goto, bypassSulfide ; Part for SULFIDES ************************************************** ; Find S-rich pixels. IF (cat[indxFeS].phsIncl EQ TRUE) THEN BEGIN if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxS,*,*] + (-1.0)*( xmaps[indxAl,*,*] + xmaps[indxSi,*,*] ) print, 'min(SmAlSi)=', MIN(r), ' max(SmAlSi)=', MAX(r) ; set values of SmACTS less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] ; for debugging: print, '...' ; histo = histogram(r, MIN=1, MAX=255) ; PLOT, histo print, '...' outFile = CritMapDir + targetfile + 'SmAlSi.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; Mark sulfides: ; Don't single out CaS + MnS + FeS + MgS sulfides in this part: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin ; Only consider pixels unknown. if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin ; Only consider the pixels that qualify as sulfide if (r[0,x,y] GE sulMinS) then begin ; Mark pixels that are sulfide, in xx map. xmaps[indxXx,x,y] = cat[indxFeS].gsval1 ; set pixel in RGB map to the RGB values for this phase. for z=0, 2 do modeRGB[z, x, y] = cat[indxFeS].rgbX[z] countSul = countSul + 1 cat[indxFeS].count = cat[indxFeS].count + 1 endif endif endfor endfor ;TV, xmaps[indxXx,*,*], /ORDER ;print, '....' ;TV, modeRGB, TRUE=1, /ORDER ;print, ' countSul = ', countSul, ' ...' outFile = CritMapDir + targetfile + 'Smode.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, modeRGB, 1 ENDIF bypassSulfide: ; goto, bypassMetal ; Part for METAL ************************************************** ; Find Fe-rich pixels. IF (cat[indxMet].phsIncl EQ TRUE) THEN BEGIN if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxFe,*,*] + (-3.0)*xmaps[indxAl,*,*] + (-3.0)*xmaps[indxCa,*,*] + (-1.0)*xmaps[indxSi,*,*] + (-1.0)*xmaps[indxMg,*,*] print, 'min(FmCMAS)=', MIN(r), ' max(FmCMAS)=', MAX(r) r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'FmCaMgSi3xAl.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; Mark metal: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] GE metMinFe) then begin xmaps[indxXx,x,y] = cat[indxMet].gsval1 ; set pixel in RGB map to the RGB values for this phase. for z=0, 2 do modeRGB[z, x, y] = cat[indxMet].rgbX[z] cat[indxMet].count = cat[indxMet].count + 1 endif endif endfor ;end loop on y endfor ;end loop on x ENDIF ; bypassMetal: ; Part for QUARTZ ************************************************** IF (cat[indxSiO2].phsIncl EQ TRUE) THEN BEGIN ; Find Si-rich pixels. if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1,maxCropElX, maxCropElY) r = xmaps[indxSi,*,*] + (-2.0) * xmaps[indxMg,*,*] + (-2.0) * xmaps[indxFe,*,*] + (-2.0) * xmaps[indxAl,*,*] + (-2.0) * xmaps[indxTi,*,*] + (-2.0) * xmaps[indxCa,*,*] ; set values of MgFeSim3Ca2AlTi less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'Sim2MgFeAl2CaTi.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; quartz ONLY: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] GE SimSiO2) then begin xmaps[indxXx,x,y] = cat[indxSiO2].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxSiO2].rgbX[z] cat[indxSiO2].count = cat[indxSiO2].count + 1 endif endif endfor endfor ENDIF ; Part for TiO2 ************************************************** IF (cat[indxTO2].phsIncl EQ TRUE) THEN BEGIN ; Find Ti-rich pixels. if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxTi,*,*] + (-2.0) * xmaps[indxMg,*,*] + (-2.0) * xmaps[indxFe,*,*] + (-2.0) * xmaps[indxAl,*,*] + (-2.0) * xmaps[indxSi,*,*] + (-2.0) * xmaps[indxCa,*,*] ; set values of Tim2MgFeAlCaSi less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'Tim2MgFeAlCaSi.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; TiO2 ONLY: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] GE TinTiO2) then begin xmaps[indxXx,x,y] = cat[indxTO2].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxTO2].rgbX[z] cat[indxTO2].count = cat[indxTO2].count + 1 endif endif endfor endfor ENDIF ; Part for CORUNDUM ************************************************** IF (cat[indxCor].phsIncl EQ TRUE) THEN BEGIN ; Find Al-rich pixels. if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxAl,*,*] + (-2.0) * xmaps[indxSi,*,*] + (-2.0) * xmaps[indxMg,*,*] + (-1.0) * xmaps[indxFe,*,*] + (-2.0) * xmaps[indxCa,*,*] ; set values of MgFeSim3Ca2AlTi less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'Alm2Si2MgFe2Ca.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; corundum ONLY: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] GE Alm2S2MF2C) then begin xmaps[indxXx,x,y] = cat[indxCor].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxCor].rgbX[z] cat[indxCor].count = cat[indxCor].count + 1 endif endif endfor endfor ENDIF ; Part for CaSiO3 WOLLASTONITE ************************************************** ; Find Ca-, Si-rich, Al-, Mg-, Fe-poor pixels. IF (cat[indxWol].phsIncl EQ TRUE) THEN BEGIN if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxCa,*,*] + xmaps[indxSi,*,*] + (-3.0) * xmaps[indxAl,*,*] + (-3.0) * xmaps[indxMg,*,*] + (-3.0) * xmaps[indxFe,*,*] ; set values of Ca3Tim3Al2Si2FePRV less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'CSm3A3M3F.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; wollastonite ONLY: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ( (xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk) ) then begin if ( r[0,x,y] GE minCSwoll ) then begin xmaps[indxXx,x,y] = cat[indxWol].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxWol].rgbX[z] cat[indxWol].count = cat[indxWol].count + 1 endif endif endfor ;end loop on y endfor ;end loop on x ENDIF ; Part for ilmenite - FeTiO3 (allow some Mg in it) ************************************************** IF (cat[indxIlm].phsIncl EQ TRUE) THEN BEGIN ; Find Mg-, Al-rich, Ca-, Si-poor pixels. if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxFe,*,*] + xmaps[indxTi,*,*] + (-2.0)*xmaps[indxSi,*,*] + (-2.0)*xmaps[indxAl,*,*] + (-2.0)*xmaps[indxCa,*,*] ; set values of MgAlm2SiCa less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'FTm2SAC.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ; TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if ( r[0,x,y] GE minFTilmenite ) AND (xmaps[indxTi,x,y] GE minTiIlm) then begin xmaps[indxXx,x,y] = cat[indxIlm].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxIlm].rgbX[z] cat[indxIlm].count = cat[indxIlm].count + 1 endif endif endfor endfor ENDIF ; Part for CaTiO3 PEROVSKITE ************************************************** ; Find Ca-, Ti-rich, Al-, Si-, Fe-poor pixels. IF (cat[indxPrv].phsIncl EQ TRUE) THEN BEGIN if (inBitDepth EQ 8) then begin r = INTARR(1, maxCropElX, maxCropElY) r = xmaps[indxCa,*,*] + 3*xmaps[indxTi,*,*] - 5*xmaps[indxAl,*,*] - xmaps[indxSi,*,*] - xmaps[indxSi,*,*] - xmaps[indxFe,*,*] - xmaps[indxFe,*,*] endif else begin r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxCa,*,*] + xmaps[indxTi,*,*] + (-1.0) * xmaps[indxAl,*,*] + (-2.0) * xmaps[indxSi,*,*] + (-2.0) * xmaps[indxFe,*,*] endelse ; set values of Ca3Tim3Al2Si2FePRV less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'Ca3Tim5Al2Si2FePRV.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; PEROVSKITE CaTiO3 ONLY: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ( (xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk) ) then begin if ( r[0,x,y] GE prvMinCaTi ) then begin xmaps[indxXx,x,y] = cat[indxPrv].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxPrv].rgbX[z] cat[indxPrv].count = cat[indxPrv].count + 1 endif endif endfor ;end loop on y endfor ;end loop on x ENDIF ; Part for HIBONITE, CA1, CA2 ************************************************** ; Choose between Aluminate minerals: ; Ca-Al-Ti only phases IF (cat[indxHib].phsIncl EQ TRUE) OR (cat[indxCA1].phsIncl EQ TRUE) OR (cat[indxCA2].phsIncl EQ TRUE) THEN BEGIN if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxCa,*,*] + (2)*xmaps[indxAl,*,*] + xmaps[indxTi,*,*] + (-3.0)*xmaps[indxSi,*,*] + (-5.0)*xmaps[indxMg,*,*] + (-3.0)*xmaps[indxFe,*,*] ; if (elms[indxS].present) then r = r + (-1.0) * xmaps[indxS,*,*] ; set values of CamS less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'CATm3S5M3F.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; Make ratio map. Add a constant to each denominator to avoid division by zero errors. (ratConstDenom = 0.01 set above) r2 = FLTARR(1, maxCropElX, maxCropElY) r2 = xmaps[indxAl,*,*] + xmaps[indxCa,*,*] - xmaps[indxMg,*,*] ;r2 = xmaps[indxAl,*,*] / ( xmaps[indxCa,*,*] + ratConstDenom ) ; r2[0,*,*] = (r[0,*,*] GT minCaCpxMelAn) * r2[0,*,*] outFile = CritMapDir + targetfile + 'AlCamMg.tif' ;outFile = CritMapDir + targetfile + 'AloverCa.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r2[0,*,*], /FLOAT ; Remove all pixels that have been identified previously: ; r2[0,*,*] = (xmaps[indxXx,*,*] EQ cat[indxUnk].gsval1) * r2[0,*,*] ; Set to zero all pixels with 1.0 < (10 * Al/Ca) < 70.0 ; r2[0,*,*] = (r2[0,*,*] GT 1.0) * r2[0,*,*] ; r2[0,*,*] = (r2[0,*,*] LT 700.0) * r2[0,*,*] for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin ; Only consider pixels unknown, AND with high Ca + Al contents. if ( (xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk) AND (r[0,x,y] GE minCATmMSF) ) then begin ; Hibonite . ;if ( (r2[0,x,y] GE AloCaminHib) AND (r2[0,x,y] LE AloCamaxHib) AND (cat[indxHib].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin if ( (r2[0,x,y] GE CAmMminHib) AND (cat[indxHib].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ; Mark pixels that are hibonite in xx map. xmaps[indxXx,x,y] = cat[indxHib].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxHib].rgbX[z] cat[indxHib].count = cat[indxHib].count + 1 endif ; CA2 CaAl2O7 (atomic Al/Ca = 1.35; wt% Al/Ca = 0.9064) if ( (r2[0,x,y] GE AloCaminCA2) AND (cat[indxCA2].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ; Mark pixels that are CA2, in xx map. xmaps[indxXx,x,y] = cat[indxCA2].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxCA2].rgbX[z] cat[indxCA2].count = cat[indxCA2].count + 1 endif ; CA1 CaAl2O4 (atomic Al/Ca = 0.67; wt% Al/Ca = 0.4532) if ( (r2[0,x,y] GE AloCaminCA1) AND (cat[indxCA1].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ; Mark pixels that are CA1, in xx map. xmaps[indxXx,x,y] = cat[indxCA1].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxCA1].rgbX[z] cat[indxCA1].count = cat[indxCA1].count + 1 endif endif ; end of conditional on unknown pixels endfor endfor r2 = 0 ; destroy r2 array. ENDIF ; Part for MgAl2O4 SPINEL ************************************************** IF (cat[indxSpn].phsIncl EQ TRUE) THEN BEGIN ; Find Mg-, Al-rich, Ca-, Si-poor pixels. if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) ;r = xmaps[indxMg,*,*] + xmaps[indxAl,*,*] + (-2.0) * xmaps[indxSi,*,*] + (-5.0) * xmaps[indxCa,*,*] ;r = (xmaps[indxCa,*,*] + xmaps[indxMg,*,*] + (-2)*xmaps[indxSi,*,*]) / (xmaps[indxAl,*,*] + ratConstDenom) ;r = (xmaps[indxAl,*,*] + xmaps[indxMg,*,*] + (-2)*xmaps[indxSi,*,*]) / (xmaps[indxCa,*,*] + ratConstDenom) r = (xmaps[indxAl,*,*] + xmaps[indxMg,*,*] + (-2)*xmaps[indxSi,*,*]) ; set values of MgAlm2SiCa less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'AMm2S.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ; TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse r2 = FLTARR(1, maxCropElX, maxCropElY) ;old r2 = xmaps[indxCa,*,*] / ( xmaps[indxAl,*,*] + xmaps[indxSi,*,*] + xmaps[indxMg,*,*] + xmaps[indxFe,*,*] + ratConstDenom ) r2 = (xmaps[indxCa,*,*] + (-1)*xmaps[indxS,*,*]) outFile = CritMapDir + targetfile + 'CamSspn.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r2[0,*,*], /FLOAT ; SPINEL ONLY: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin ;if ( r[0,x,y] GE MgAlm2Si5CaminSpn ) then begin ;if ( r[0,x,y] GE CM2SoverAminSpn ) then begin if ( r[0,x,y] GE AMm2SminSpn ) AND (r[0,x,y] LE AMm2SmaxSpn) AND (r2[0,x,y] LE maxCaspn) then begin xmaps[indxXx,x,y] = cat[indxSpn].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxSpn].rgbX[z] cat[indxSpn].count = cat[indxSpn].count + 1 endif endif endfor endfor r2 = 0 ENDIF ; Part for sphene-CaTiSiO5 ************************************************** IF (cat[indxTtn].phsIncl EQ TRUE) THEN BEGIN ; Find Mg-, Al-rich, Ca-, Si-poor pixels. if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxCa,*,*] + xmaps[indxTi,*,*] + xmaps[indxSi,*,*] + (-5.0) * xmaps[indxAl,*,*] + (-5.0) * xmaps[indxMg,*,*] + (-5.0) * xmaps[indxFe,*,*] ; set values of MgAlm2SiCa less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'CTSm5AMF.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ; TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if ( r[0,x,y] GE minCTSsphene ) AND (xmaps[indxTi,x,y] GE minTisphene) then begin xmaps[indxXx,x,y] = cat[indxTtn].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxTtn].rgbX[z] cat[indxTtn].count = cat[indxTtn].count + 1 endif endif endfor endfor ENDIF ; Part for CARBONATE ************************************************** IGNORES Dolomite !!! ; Find Ca-rich pixels w/ no other elements. IF (cat[indxCc].phsIncl EQ TRUE) THEN BEGIN if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = (xmaps[indxCa,*,*] + (-2.0)*(xmaps[indxAl,*,*] + xmaps[indxSi,*,*] + xmaps[indxMg,*,*] + xmaps[indxFe,*,*])) r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'Cm2ASMF.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; Mark carbonate: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin ; Only consider pixels unknown. if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] GE ccMinCa) then begin ; Mark pixels that are carbonate, in xx map. xmaps[indxXx,x,y] = cat[indxCc].gsval1 ; set pixel in RGB map to the RGB values for this phase. for z=0, 2 do modeRGB[z, x, y] = cat[indxCc].rgbX[z] cat[indxCc].count = cat[indxCc].count + 1 endif endif endfor endfor ENDIF ; Part for Ca and Al rich pyroxene (Ca-Tschermak's or fassaite)************************************************** IF (cat[indxCaTs].phsIncl EQ TRUE) THEN BEGIN ; Find Ca-, Al-rich pixels. if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1,maxCropElX, maxCropElY) r = (xmaps[indxAl,*,*]+ xmaps[indxCa,*,*] ) / (xmaps[indxSi,*,*] + ratConstDenom) ; set values less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'ACoverS.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ; TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if ( r[0,x,y] LE ACoverSmaxCaTs ) AND ( r[0,x,y] GE ACoverSminCaTs ) then begin xmaps[indxXx,x,y] = cat[indxCaTs].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxCaTs].rgbX[z] cat[indxCaTs].count = cat[indxCaTs].count + 1 endif endif endfor endfor ENDIF ; Part for FeAl2O4 SPINEL AKA Hercynite************************************************** IF (cat[indxHrc].phsIncl EQ TRUE) THEN BEGIN ; Find Fe-, Al-rich, Ca-, Si-poor pixels. if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = (xmaps[indxAl,*,*] + xmaps[indxFe,*,*] + (-2)*xmaps[indxSi,*,*]) / (xmaps[indxCa,*,*] + ratConstDenom) ; set values of MgAlm2SiCa less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'AFm2SoverC.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ; TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin ;if ( r[0,x,y] GE MgAlm2Si5CaminSpn ) then begin ;if ( r[0,x,y] GE CM2SoverAminSpn ) then begin if ( r[0,x,y] GE AFm2SoverCminHrc ) then begin xmaps[indxXx,x,y] = cat[indxHrc].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxHrc].rgbX[z] cat[indxHrc].count = cat[indxHrc].count + 1 endif endif endfor endfor ENDIF ; Part for OLIVINE and OPX ************************************************** IF (cat[indxOlv].phsIncl EQ TRUE) OR (cat[indxOpx].phsIncl EQ TRUE) THEN BEGIN ; Find Mg-Fe-Si rich pixels if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxMg,*,*] + xmaps[indxFe,*,*] + xmaps[indxSi,*,*] + (-3.0) * xmaps[indxCa,*,*] + (-2.0) * xmaps[indxAl,*,*] + (-1.0) * xmaps[indxTi,*,*] ; set values of MgFeSim3Ca2AlTi less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'MgFeSim3Ca2AlTi.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse r2 = FLTARR(1, maxCropElX, maxCropElY) ; Make ratio map. Add a constant to each denominator to avoid division by zero errors. (ratConstDenom = 0.01 set above) r2 = ( xmaps[indxMg,*,*] + xmaps[indxFe,*,*] ) / ( xmaps[indxSi,*,*] + ratConstDenom ) r2[0,*,*] = (r[0,*,*] GE MgFeSim3Ca2AlTiminOLOPX) * r2[0,*,*] outFile = CritMapDir + targetfile + 'MgFeoverSi.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r2[0,*,*], /FLOAT for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin ; Find Mg-Fe-Si rich pixels and then use ratio (Mg+Fe)/Si to tell OL and OPX apart. if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] GE MgFeSim3Ca2AlTiminOLOPX) AND (r2[0,x,y] GE MgFevsSiPX) then begin if (r2[0,x,y] GE MgFevsSiOLmin) AND (r2[0,x,y] LE MgFevsSiOLmax) AND (cat[indxOlv].phsIncl EQ TRUE) then begin xmaps[indxXx,x,y] = cat[indxOlv].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxOlv].rgbX[z] cat[indxOlv].count = cat[indxOlv].count + 1 endif else begin if (cat[indxOpx].phsIncl EQ TRUE) then begin xmaps[indxXx,x,y] = cat[indxOpx].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxOpx].rgbX[z] cat[indxOpx].count = cat[indxOpx].count + 1 endif endelse endif endif endfor ; end loop on y. endfor ; end loop on x r2 = 0 ; destroy r2 array. ENDIF goto, DSEplag ; ECP Part for plagioclase (anorthite) and other feldspar. IF (cat[indxPlg].phsIncl EQ TRUE) THEN BEGIN if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = (5)*xmaps[indxSi,*,*] +(-1)*xmaps[indxAl,*,*] +(-1)*xmaps[indxCa,*,*] +(-1)*xmaps[indxMg,*,*] ; set values of 5SmACM less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + '5SmACM.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ; TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse r2 = FLTARR(1, maxCropElX, maxCropElY) ; Make ratio map. Add a constant to each denominator to avoid division by zero errors. (ratConstDenom = 0.01 set above) r2 = (xmaps[indxSi,*,*] + xmaps[indxAl,*,*])/ (2*xmaps[indxCa,*,*] + xmaps[indxMg,*,*] + ratConstDenom ) r2[0,*,*] = (r[0,*,*] GT minSiplg) * r2[0,*,*] outFile = CritMapDir + targetfile + 'SAover2CM.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r2[0,*,*], /FLOAT for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] GE minSiplg)and (r2[0,x,y] GE minplag ) then begin xmaps[indxXx,x,y] = cat[indxPlg].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxPlg].rgbX[z] cat[indxPlg].count = cat[indxPlg].count + 1 endif endif endfor endfor r2 = 0 ; destroy r2 array. ENDIF ;goto, bypassDSEplag DSEplag: IF (cat[indxPlg].phsIncl EQ TRUE) THEN BEGIN if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxSi,*,*] + 4*xmaps[indxAl,*,*] + xmaps[indxCa,*,*] + (-10.0*(xmaps[indxMg,*,*]) + xmaps[indxFe,*,*] + (2*xmaps[indxTi,*,*])) if (elms[indxNa].present) then r = r + xmaps[indxNa,*,*] ; set values of SCAm10MFT less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'SCAm10MFT.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ; TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse r2 = FLTARR(1, maxCropElX, maxCropElY) ; Make ratio map. Add a constant to each denominator to avoid division by zero errors. (ratConstDenom = 0.01 set above) r2 = xmaps[indxSi,*,*]/ ( 5*xmaps[indxAl,*,*] + ratConstDenom ) outFile = CritMapDir + targetfile + 'SioverAl.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r2[0,*,*], /FLOAT for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin ;if (r[0,x,y] GE minSCAplg) then begin if (r[0,x,y] GE minSCAplg) then begin if (r2[0,x,y] GE minSoAplag ) and (r2[0,x,y] LE maxSoAplag ) then begin xmaps[indxXx,x,y] = cat[indxPlg].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxPlg].rgbX[z] cat[indxPlg].count = cat[indxPlg].count + 1 endif else begin if (elms[indxNa].present) then begin if (xmaps[indxNa,x,y] GE minNaAlbite) AND (r2[0,x,y] LE maxSoAalb) then begin for z=0, 2 do modeRGB[z, x, y] = cat[indxAb].rgbX[z] cat[indxAb].count = cat[indxAb].count + 1 endif endif endelse endif endif endfor endfor r2 = 0 ; destroy r2 array. ENDIF bypassDSEplag: ; Find Ca-rich phases: CPX MELILITE : IF (cat[indxCpx].phsIncl EQ TRUE) OR (cat[indxMel].phsIncl EQ TRUE) OR (cat[indxGls].phsIncl EQ TRUE) OR (cat[indxPlg].phsIncl EQ TRUE) THEN BEGIN if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxCa,*,*] ;+ (-1.0) * xmaps[indxFe,*,*] if (elms[indxS].present) then r = r + (-1.0) * xmaps[indxS,*,*] ;+ (-1.0) *xmaps[indxMg,*,*] + (-1.0) *xmaps[indxSi,*,*] + (-2.0) *xmaps[indxFe,*,*] + (-1.0) *xmaps[indxTi,*,*] ; set values of CamS less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'CamS.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse ; Make ratio map. Add a constant to each denominator to avoid division by zero errors. (ratConstDenom = 0.01 set above) r2 = FLTARR(1, maxCropElX, maxCropElY) ;old r2 = xmaps[indxCa,*,*] / ( xmaps[indxAl,*,*] + xmaps[indxSi,*,*] + xmaps[indxMg,*,*] + xmaps[indxFe,*,*] + ratConstDenom ) r2 = (xmaps[indxSi,*,*]) / (xmaps[indxCa,*,*] + ratConstDenom ) ;r2[0,*,*] = (r[0,*,*] GT minCaCpxMelAn) * r2[0,*,*] outFile = CritMapDir + targetfile + 'SoC.tif' ;r2 = (2*xmaps[indxCa,*,*]) / (xmaps[indxSi,*,*] + xmaps[indxMg,*,*] + ratConstDenom ) ;r2[0,*,*] = (r[0,*,*] GT minCaCpxMelAn) * r2[0,*,*] ;outFile = CritMapDir + targetfile + '2CoverMS.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r2[0,*,*], /FLOAT ;r3 = FLTARR(1, maxCropElX, maxCropElY) ;r3 = (xmaps[indxSi,*,*]) + (xmaps[indxAl,*,*]) + (2.0) * (xmaps[indxMg,*,*]) + (-2.0) * (xmaps[indxCa,*,*]) ;outFile = CritMapDir + targetfile + 'SiAl2Mgm2Ca.tif' r3 = FLTARR(1, maxCropElX, maxCropElY) r3 = (xmaps[indxSi,*,*]) + (xmaps[indxAl,*,*]) - (xmaps[indxCa,*,*]) outFile = CritMapDir + targetfile + 'SiAlmCa.tif' ;r3 = FLTARR(1, maxCropElX, maxCropElY) ;r3 = (xmaps[indxMg,*,*]) + (xmaps[indxSi,*,*]) + (xmaps[indxCa,*,*]) ;outFile = CritMapDir + targetfile + 'MgSiCa.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r3[0,*,*], /FLOAT for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin ; Only consider pixels unknown, AND with high Ca ratio. if ( (xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk) AND (r[0,x,y] GE minCaCpxMelAn) ) then begin ; Melilite (Mel) has C/MAS = 0.667 atoms. ;if ( (r2[0,x,y] LE SoCmaxmel) AND (r2[0,x,y] GE SoCminmel) AND (r3[0,x,y] LE SAMCmaxmel) AND (cat[indxMel].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin if ( (r2[0,x,y] LE SoCmaxmel) AND (r2[0,x,y] GE SoCminmel) AND (cat[indxMel].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ;if ( (r2[0,x,y] GE CoMSminmel) AND (cat[indxMel].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ;if ( (r2[0,x,y] GE CaMinCoMASFmel) AND (cat[indxMel].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ; Mark pixels that are melilite in xx map. xmaps[indxXx,x,y] = cat[indxMel].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxMel].rgbX[z] cat[indxMel].count = cat[indxMel].count + 1 endif ; CPX (fassaite) has C/MAS = 1.0 atoms. ;if ( (r2[0,x,y] GE SoCmincpx) AND (r2[0,x,y] LE SoCmaxcpx) AND (r3[0,x,y] GE SAMCmincpx) AND (r3[0,x,y] LE SAMCmaxcpx) AND (cat[indxCpx].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin if ( (r2[0,x,y] GE SoCmincpx) AND (r2[0,x,y] LE SoCmaxcpx) AND (r3[0,x,y] LE SACmaxcpx) AND (cat[indxCpx].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ;if ( (r2[0,x,y] GE SoCmincpx) AND (r2[0,x,y] LE SoCmaxcpx) AND (r3[0,x,y] GE MSCmincpx) AND (cat[indxCpx].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ;if ( (r2[0,x,y] GE CoMSmincpx) AND (cat[indxCpx].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ;if ( (r2[0,x,y] GE CaMinCoMASFcpx) AND (cat[indxCpx].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin ; Mark pixels that are Cpx, in xx map. xmaps[indxXx,x,y] = cat[indxCpx].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxCpx].rgbX[z] cat[indxCpx].count = cat[indxCpx].count + 1 endif endif if ( (xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk) AND (r[0,x,y] GE minCameso) ) then begin if ( (r2[0,x,y] GE SoCminmeso) AND (r2[0,x,y] LE SoCmaxmeso) AND (cat[indxGls].phsIncl EQ TRUE) AND (xmaps[indxXx,x,y] EQ gsUnk) ) then begin xmaps[indxXx,x,y] = cat[indxGls].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxGls].rgbX[z] cat[indxGls].count = cat[indxGls].count + 1 endif endif ; end of conditional on unknown pixels endfor endfor r2 = 0 ; destroy r2 array. r3 = 0 ENDIF goto, bypassNephAbSod ; Na-, Al-silicates Ca-free: Alteration sodalite Na8(AlSiO4)6Cl2 or nepheline NaAlSiO4 or albite NaAlSi3O8: if (inBitDepth EQ 8) then r = INTARR(1, maxCropElX, maxCropElY) else r = FLTARR(1, maxCropElX, maxCropElY) r = xmaps[indxAl,*,*] + xmaps[indxSi,*,*] + (-3.0) * xmaps[indxCa,*,*] + (-3.0) * xmaps[indxMg,*,*] + (-3.0) * xmaps[indxFe,*,*] ; set values of Alm3Ca3Mg3FeminAlt less than 0, to zero value. r[0,*,*] = (r[0,*,*] GT 0.0) * r[0,*,*] outFile = CritMapDir + targetfile + 'Alm3Ca3Mg3Fe.tif' if (inBitDepth EQ 8) then begin ; Scale r to 0-255: r[0,*,*] = BYTSCL(r[0,*,*]) ;TV, r[0,*,*], /ORDER ;print, '...' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], 1 endif else begin if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, r[0,*,*], /FLOAT endelse for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if ((xmaps[indxXx,x,y] EQ gsUnk) AND (xmaps[indxXx,x,y] NE gsMsk)) then begin if (r[0,x,y] GE altMinNa) then begin xmaps[indxXx,x,y] = cat[indxAlt].gsval1 for z=0, 2 do modeRGB[z, x, y] = cat[indxAlt].rgbX[z] cat[indxAlt].count = cat[indxAlt].count + 1 endif endif endfor endfor bypassNephAbSod: ; end of phase decisions. ; Now calculate the number of unknown pixels in the masked and cropped area: for x=0, (maxCropElX-1) do begin for y=0, (maxCropElY-1) do begin if (xmaps[indxXx,x,y] EQ cat[indxUnk].gsval1) then cat[indxUnk].count = cat[indxUnk].count + 1 endfor endfor ; Calculate the number of masked pixels: countObj = 0L countTot = 0L ; add 3july2016 DSE for p=0, (sizePhaseArray-1) do begin if (p NE indxMsk) then countObj = countObj + cat[p].count endfor countTot = LONG(maxCropElX)*LONG(maxCropElY) ; add 3july2016 DSE. countTot must be a LONG integer! cat[indxMsk].count = countTot - countObj ENDIF ; end of conditional on MapCrit EQ 2 ; Mode has been determined ======================================================================= ; Now calculate mean intensity/pixel and its std-dev for each element in each phase: ; Loop through all included phases. For each phase, run through mode map. Sum element counts in each phase. ; Loop over all pixels. ; 6-July-2016 (DSE) Bug correction! x was used twice as variable in this segment! ; cat[].elemSums[] are set to 0.0 above. for x=0, (maxCropElX-1) do begin ; loop through all x-y in the maps for y=0, (maxCropElY-1) do begin for p=0, (sizePhaseArray-1) do begin ; loop through all phases p, if cat[p].phsIncl then begin ; do NOT include all phases ; Find what phase this pixel is, then sum map values for each element, and break out of phase loop. if (xmaps[indxXx,x,y] EQ cat[p].gsval1) then begin for z=0, (numImages-1) do begin cat[p].elemSums[z] = cat[p].elemSums[z] + DOUBLE(xmaps[z,x,y]) ; make DBL 12july2016 DSE elemAvgVal[z,5] = elemAvgVal[z,5] + DOUBLE(xmaps[z,x,y]) endfor break endif endif endfor endfor endfor ; set value for curSampleN from cat[].elemSums[] for p=0, (sizePhaseArray-1) do begin for z=1, (numImages-1) do begin allelmSum[p, curSampleN-1, z] = cat[p].elemSums[z] ; The next item causes a 'divide by zero' error if (cat[p].count GT 0) then cat[p].elemMeans[z] = cat[p].elemSums[z] / cat[p].count allelmMean[p, curSampleN-1, z] = cat[p].elemMeans[z] endfor endfor ; ================================================================== ; NOW OUTPUT RESULT ; ================================================================== ; Print TIF of the mode map: outFile = CritMapDir + targetfile + '_RGBmode.tif' WRITE_TIFF, outFile, modeRGB, 1 outFile = CritMapDir + targetfile + '_GSmode.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, xmaps[indxXx,*,*], /FLOAT bypassModeMapping: ; Print TIF of the SOURCE elements, masked, as RGB maps: print, ' BEGIN: Construction of 3-color mosaics.' RGBimg = MAKE_ARRAY(3, maxCropElX, maxCropElY, /BYTE) ; Autolevel all the 32-bit maps (we are done with them at this point). maxCutoff = 0.995 minCutoff = 0.995 tempMap = BYTARR(numImages, maxCropElX, maxCropElY) for n=0, (numImages-1) do begin ; make 8-bit arrays out of the xmaps tempMap[n, 0:(maxCropElX-1), 0:(maxCropElY-1)] = BYTSCL(xmaps[n,0:(maxCropElX-1), 0:(maxCropElY-1)]) ; tempMap[n, 0:(maxCropElX-1), 0:(maxCropElY-1)] = BYTSCL(tempMap, TOP=255) ) ; calls the procedure stored at path: E;\a\IDL-libs ; CALL SUBPROGRAM auto_Level_v1, elMap, minCutoff, maxCutoff < arguments if (inBitDepth EQ 32) then auto_Level_v2, tempMap[n, *, *], minCutoff, maxCutoff endfor ; Mg-Ca-Al if (elms[indxMg].present AND elms[indxCa].present AND elms[indxAl].present) then begin RGBimg[0, *, *] = tempMap[indxMg, *, *] RGBimg[1, *, *] = tempMap[indxCa, *, *] RGBimg[2, *, *] = tempMap[indxAl, *, *] outFile = CritMapDir + targetfile + '_MgCaAl.tiff' print, 'Writing full Mg-Ca-Al RGB mosaic # to ', outFile,'.' WRITE_TIFF, outFile, RGBimg, ORIENTATION=1 endif ; Ti-Ca-Al if (elms[indxTi].present AND elms[indxCa].present AND elms[indxAl].present) then begin RGBimg[0, *, *] = tempMap[indxTi, *, *] RGBimg[1, *, *] = tempMap[indxCa, *, *] RGBimg[2, *, *] = tempMap[indxAl, *, *] outFile = CritMapDir + targetfile + '_TiCaAl.tiff' print, 'Writing full Ti-Ca-Al RGB mosaic # to ', outFile,'.' WRITE_TIFF, outFile, RGBimg, ORIENTATION=1 endif ; Si-Ca-Fe if (elms[indxSi].present AND elms[indxCa].present AND elms[indxFe].present) then begin RGBimg[0, *, *] = tempMap[indxSi, *, *] RGBimg[1, *, *] = tempMap[indxCa, *, *] RGBimg[2, *, *] = tempMap[indxFe, *, *] outFile = CritMapDir + targetfile + '_SiCaFe.tiff' print, 'Writing full Si-Ca-Fe RGB mosaic # to ', outFile,'.' WRITE_TIFF, outFile, RGBimg, ORIENTATION=1 endif ; Print a tiff file with bars showing colors of all the items (the legend) ;If (numSamples EQ filenum-1) then begin width = 40 ;replace 22jun2016 width = FIX( ((maxCropElY-1)/sizePhaseArray)-1 ) heightInc = 20 height = 0 i=0 for x=0, sizePhaseArray-1 do begin if (cat[x].phsIncl EQ 1) then height = height+heightInc if (cat[x].phsIncl EQ 1) then i = i+1 endfor legendRGB = BYTARR(3, width, height+1) print, 'ledgend width=', width, ' height= ', height, ' numPhases=', i, ' sizePhaseArray = ', sizePhaseArray ; loop thru all the phases that are present. Make a bar of fixed width for each color. i = 0 for x=0, sizePhaseArray-1 do begin if (cat[x].phsIncl EQ 1) then begin start = i * heightInc + 2 i = i+1 for y=start, (start+heightInc-4) do begin for z=0, 2 do begin legendRGB[z, *, y] = cat[x].rgbX[z] endfor endfor endif endfor outFile = ModeDir + targetfile + '_legend.tif' if (writeIntermedTifYes1no2 EQ 1) then WRITE_TIFF, outFile, legendRGB, 1 print, 'Phases present (LEGEND)' for x=0, sizePhaseArray-1 do if (cat[x].phsIncl EQ 1) then print, cat[x].label ; Print information to LOG file: if (countSul EQ 0) then countSul = cat[indxFeS].count + cat[indxCaS].count + cat[indxMnS].count + cat[indxMgS].count ; cat[indxUnk].count = 0 < this is already calculated above. unkCheck = 0 for x=0, sizePhaseArray-1 do begin if (x NE indxUnk) AND (x NE indxMsk) then unkCheck = unkCheck + cat[x].count ;Could add AND (x NE indxHole) endfor unkCheck = countObj - unkCheck if (unkCheck NE cat[indxUnk].count) then print, 'WARNING: Unknown value NOT correct' ;ENDIF ;printf, loglun, '==================================================================' ; Print avg & stdev of pixel values for each phase included in the mode. ;printf, loglun, '==================================================================' ;printf, loglun, 'date:', STRING(9B), SYSTIME() ; Print numbers of pixel, etc. for each phase included in the mode. print, 'date: ', STRING(9B),SYSTIME() print, 'sample: ', STRING(9B),basefile print, 'size: ', STRING(9B), maxCropElX,STRING(9B),maxCropElY, STRING(9B), ' pixels', STRING(9B), maxCropElX*maxCropElY, STRING(9B), '