source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/SANS/Marquee.ipf @ 575

Last change on this file since 575 was 570, checked in by srkline, 13 years ago

Change (1):
In preparation for release, updated pragma IgorVersion?=6.1 in all procedures

Change (2):
As a side benefit of requiring 6.1, we can use the MultiThread? keyword to thread any model function we like. The speed benefit is only noticeable on functions that require at least one integration and at least 100 points (resolution smearing is NOT threaded, too many threadSafe issues, too little benefit). I have chosen to use the MultiThread? only on the XOP assignment. In the Igor code there are too many functions that are not explicitly declared threadsafe, making for a mess.

File size: 15.5 KB
Line 
1#pragma rtGlobals=1             // Use modern global access method.
2#pragma version=5.0
3#pragma IgorVersion=6.1
4
5////////////
6// vers 1.21 3 may 06 total transmission incorporated (BSG)
7//
8//**************************
9// Vers 1.2 091901
10//
11// marquee functions are used to:
12//
13// locate the beamcenter
14// set the (x,y) box coordinates for the sum for transmission calculation
15// read out coordinates
16// do a "box sum" of the same box over a range of files
17// do a 2D Gaussian fit over a selected range
18// do a save of the current image (with colorBar) - as a graphics image
19//
20//***************************
21
22//sums the data counts in the box specified by (x1,y1) to (x2,y2)
23//assuming that x1<x2, and y1<y2
24//the x,y values must also be in axis coordinates[0,127] NOT (1,128) detector coords.
25//
26// accepts arbitrary detector coordinates. calling function is responsible for
27// keeping selection in bounds
28Function SumCountsInBox(x1,x2,y1,y2,type)
29        Variable x1,x2,y1,y2
30        String type
31       
32        Variable counts = 0,ii,jj
33       
34        String dest =  "root:Packages:NIST:"+type
35       
36        //check for logscale data, but don't change the data
37        NVAR gIsLogScale = $(dest + ":gIsLogScale")
38        if (gIsLogScale)
39                wave w=$(dest + ":linear_data")
40        else
41                wave w=$(dest + ":data")
42        endif
43       
44        ii=x1
45        jj=y1
46        do
47                do
48                        counts += w[ii][jj]
49                        jj+=1
50                while(jj<=y2)
51                jj=y1
52                ii+=1
53        while(ii<=x2)
54       
55        Return (counts)
56End
57
58
59//from a marquee selection:
60//calculates the sum of counts in the box, and records this count value in
61//globals for both patch and Trans routines, then records the box coordinates
62//and the count value for the box in the header of the (currently displayed)
63//empty beam file for use later when calculating transmissions of samples
64//these values are written to unused (analysis) fields in the data header
65//4 integers and one real value are written
66//
67//re-written to work with Transmission panel as well as PatchPanel
68//which now bothe work from the same folder (:Patch)
69//
70Function SetXYBoxCoords() :  GraphMarquee
71
72        GetMarquee left,bottom
73        if(V_flag == 0)
74                Abort "There is no Marquee"
75        Endif
76        SVAR dispType=root:myGlobals:gDataDisplayType
77        if(cmpstr(dispType,"SAM")!=0)
78                DoAlert 0, "You can only use SetXYBox on SAM data files"
79                return(1)
80        endif
81        //printf "marquee left in bottom axis terms: %g\r",round(V_left)
82        //printf "marquee right in bottom axis terms: %g\r",round(V_right)
83        //printf "marquee top in left axis terms: %g\r",round(V_top)
84        //printf "marquee bottom in left axis terms: %g\r",round(V_bottom)
85        Variable x1,x2,y1,y2
86        x1 = round(V_left)
87        x2 = round(V_right)
88        y1 = round(V_bottom)
89        y2 = round(V_top)
90       
91        KeepSelectionInBounds(x1,x2,y1,y2)
92       
93        //check to make sure that Patch and Trans data folders exist for writing of global variables
94        If( ! (DataFolderExists("root:myGlobals:Patch"))  )
95                Execute "InitializePatchPanel()"
96        Endif
97        //check to make sure that Patch and Trans data folders exist for writing of global variables
98        If( ! (DataFolderExists("root:myGlobals:TransHeaderInfo"))  )
99                Execute "InitializeTransPanel()"
100        Endif
101       
102        //write string as keyword-packed string, to use IGOR parsing functions
103        String msgStr = "X1="+num2str(x1)+";"
104        msgStr += "X2="+num2str(x2)+";"
105        msgStr += "Y1="+num2str(y1)+";"
106        msgStr += "Y2="+num2str(y2)+";"
107        String/G root:myGlobals:Patch:gPS3 = msgStr
108        String/G root:myGlobals:Patch:gEmpBox = msgStr
109        //changing this global wil update the display variable on the TransPanel
110        String/G root:myGlobals:TransHeaderInfo:gBox = msgStr
111       
112        //sum the counts in the patch - working on the SAM data, to be sure that it's normalized
113        //to the same monitor counts and corrected for detector deadtime
114        String type = "SAM"
115        Variable counts
116        counts = SumCountsInBox(x1,x2,y1,y2,type)
117        Print " marquee counts =",counts
118        //Set the global gTransCts
119        Variable/G root:myGlobals:Patch:gTransCts = counts
120       
121        //now change the extra variables in the empty beam file
122        //get the filename from the SAM folder (there will only be one file)
123        SVAR partialName = root:Packages:NIST:SAM:FileList
124        //construct valid filename, then prepend path
125        String tempName = FindValidFilename(partialName)
126        Print "in marquee",partialName
127        //Print tempName
128        if(cmpstr(tempName,"")==0)
129                //file not found, get out
130                Abort "file not found, marquee"
131        Endif
132        //name is ok, prepend path to tempName for read routine
133        PathInfo catPathName
134        if (V_flag == 0)
135                //path does not exist - no folder selected
136                Abort "no path selected"
137        else
138                String filename = S_path + tempName
139        endif
140       
141        if(cmpstr(filename,"no file selected")==0)
142                Abort "no file selected"
143        Endif
144       
145        WriteXYBoxToHeader(filename,x1,x2,y1,y2)
146       
147        Print counts, " counts in XY box"
148        WriteBoxCountsToHeader(filename,counts)
149       
150End
151
152//finds the beam center (the centroid) of the selected region
153//and simply prints out the results to the history window
154//values are printed out in detector coordinates, not IGOR coords.
155//
156Function FindBeamCenter() :  GraphMarquee
157
158        //get the current displayed data (so the correct folder is used)
159        SVAR cur_folder=root:myGlobals:gDataDisplayType
160        String dest = "root:Packages:NIST:" + cur_folder
161       
162        Variable xzsum,yzsum,zsum,xctr,yctr
163        Variable left,right,bottom,top,ii,jj,counts
164       
165        // data wave is hard-wired in as the displayed data
166        NVAR dataIsLog=$(dest + ":gIsLogScale")         //check for log-scaling in current data folder
167        if (dataIsLog)
168                wave data=$(dest + ":linear_data")
169        else
170                wave data=$(dest + ":data")
171        endif
172       
173        GetMarquee left,bottom
174        if(V_flag == 0)
175                Print "There is no Marquee"
176        else
177                left = round(V_left)            //get integer values for selection limits
178                right = round(V_right)
179                top = round(V_top)
180                bottom = round(V_bottom)
181               
182                KeepSelectionInBounds(left,right,bottom,top)
183               
184                // selection valid now, calculate beamcenter
185                xzsum = 0
186                yzsum = 0
187                zsum = 0
188                // count over rectangular selection, doing each row, L-R, bottom to top
189                ii = bottom -1
190                do
191                        ii +=1
192                        jj = left-1
193                        do
194                                jj += 1
195                                counts = data[jj][ii]
196                                xzsum += jj*counts
197                                yzsum += ii*counts
198                                zsum += counts
199                        while(jj<right)
200                while(ii<top)
201               
202                xctr = xzsum/zsum
203                yctr = yzsum/zsum
204               
205                // add 1 to each to get to detector coordinates (1,128)
206                // rather than the data array which is [0,127]
207                xctr+=1
208                yctr+=1
209               
210                Print "X-center (in detector coordinates) = ",xctr
211                Print "Y-center (in detector coordinates) = ",yctr
212        endif
213       
214        //back to root folder (redundant)
215        SetDataFolder root:
216       
217End
218
219//still need to error check - out-of-bounds...waves exist.
220// allows a 2D Gaussian fit to a selected region of data in a SANS_Data window
221//puts up a new graph with the fitted contour
222Function Do_2D_Gaussian_Fit() :  GraphMarquee
223        String topWin=WinName(0,1)              //top *graph* window
224        //exit nicely if not in the Data display window
225        if(cmpstr(topWin,"SANS_Data") != 0)
226                DoAlert 0,"2D Gaussian fitting is only available from the Data Display Window"
227                return(1)
228        Endif
229       
230        GetMarquee/K left,bottom
231        Variable x1,x2,y1,y2,qxlo,qxhi,qylo,qyhi
232        if(V_flag == 0)
233                Print "There is no Marquee"
234        else
235                String junk="",df=""
236               
237                //**hard-wired info about the x-y q-scales
238                qxlo = DimOffset(root:myGlobals:q_x_axis,0)
239                qxhi = DimDelta(root:myGlobals:q_x_axis,0) + qxlo
240//              Print "qxlo,qxhi = ",qxlo,qxhi
241                Wave w=$"root:myGlobals:q_y_axis"
242                qylo=w[0]
243                qyhi=w[1]
244//              print "qylo,qyhi = ",qylo,qyhi
245               
246                junk=ImageInfo("SANS_Data","data",0)
247                df=StringByKey("ZWAVEDF", junk,":",";")
248//              print df
249                Duplicate/O $(df+"data") data,data_err
250                data_err=sqrt(data)             //for weighting
251               
252// comment out the SetScale lines if you want the result in terms of pixels as a way of
253// measuring the beam center. Note that you need to ADD ONE to fitted x0 and y0 to get detector
254// coordinates rather than the zero-indexed array. 2D fitting does have the benefit of
255// reporting error bars on the xy (if you believe that 2D gaussian is correct)         
256                SetScale/I x qxlo,qxhi,"",data
257                SetScale/I y qylo,qyhi,"",data
258               
259                Display /W=(10,50,361,351) /K=1
260                AppendImage data
261                ModifyImage data ctab= {*,*,Grays,0}
262                ModifyGraph width={Plan,1,bottom,left}
263                ModifyGraph mirror=2
264                ModifyGraph lowTrip=1e-04
265                ModifyImage data cindex=$"root:myGlobals:NIHColors"
266                SVAR angst = root:myGlobals:gAngstStr
267                Label bottom "Qx ("+angst+"\\S-1\\M)"
268                Label left "Qy ("+angst+"\\S-1\\M)"
269
270                //keep selection in-bounds
271                x1=V_left
272                x2=V_right
273                y1=V_bottom
274                y2=V_top
275                KeepSelectionInBounds(x1,x2,y1,y2)
276
277                //cross correlation coefficent (K6) must be between 0 and 1, need constraints
278                Make/O/T/N=2 temp_constr
279                temp_constr = {"K6>0","K6<1"}
280               
281                CurveFit/N Gauss2D data[x1,x2][y1,y2] /I=1 /W=data_err /D /R /A=0 /C=temp_constr
282               
283                Killwaves/Z temp_constr
284        endif
285End
286
287// to save the image, simply invoke the IGOR menu item for saving graphics
288//
289Function SaveSANSGraphic() : GraphMarquee
290       
291        NVAR isDemoVersion=root:myGlobals:isDemoVersion
292        if(isDemoVersion==1)
293                //      comment out in DEMO_MODIFIED version, and show the alert
294                DoAlert 0,"This operation is not available in the Demo version of IGOR"
295        else
296                DoIGORMenu "File","Save Graphics"
297        endif
298End
299
300//does a sum over each of the files in the list over the specified range
301// x,y are assumed to already be in-bounds of the data array
302// output is dumped to the command window
303//
304Function DoBoxSum(fileStr,x1,x2,y1,y2,type)
305        String fileStr
306        Variable x1,x2,y1,y2
307        String type
308       
309        //parse the list of file numbers
310        String fileList="",item="",pathStr="",fullPath=""
311        Variable ii,num,err,cts
312       
313        PathInfo catPathName
314        If(V_Flag==0)
315                Abort "no path selected"
316        Endif
317        pathStr = S_Path
318       
319        fileList=ParseRunNumberList(fileStr)
320        num=ItemsInList(fileList,",")
321       
322        //loop over the list
323        //add each file to SAM (to normalize to monitor counts)
324        //sum over the box
325        //print the results
326        Make/O/N=(num) FileID,BoxCounts
327        Print "Results are stored in root:FileID and root:BoxCounts waves"
328        for(ii=0;ii<num;ii+=1)
329                item=StringFromList(ii,fileList,",")
330                FileID[ii] = GetRunNumFromFile(item)            //do this here, since the list is now valid
331                fullPath = pathStr+item
332                ReadHeaderAndData(fullPath)
333//              String/G root:myGlobals:gDataDisplayType="RAW"
334//              fRawWindowHook()
335                if(cmpstr(type,"SAM")==0)
336                        err = Raw_to_work("SAM")
337                endif
338                String/G root:myGlobals:gDataDisplayType=type
339                fRawWindowHook()
340                cts=SumCountsInBox(x1,x2,y1,y2,type)
341                BoxCounts[ii]=cts
342                Print item+" counts = ",cts
343        endfor
344       
345        DoBoxGraph(FileID,BoxCounts)
346       
347        return(0)
348End
349
350Function DoBoxGraph(FileID,BoxCounts)
351        Wave FileID,BoxCounts
352       
353        Sort FileID BoxCounts,FileID            //sort the waves, in case the run numbers were entered out of numerical order
354       
355        Display /W=(5,44,383,306) BoxCounts vs FileID
356        ModifyGraph mode=4
357        ModifyGraph marker=8
358        ModifyGraph grid=2
359        ModifyGraph mirror=2
360        Label left "Counts (per 10^8 monitor counts)"
361        Label bottom "Run Number"
362       
363        return(0)
364End
365
366//
367// promts the user for a range of file numbers to perform the sum over
368// list must be comma delimited numbers (or dashes) just as in the BuildProtocol panel
369// the (x,y) range is already selected from the marquee
370//
371Function BoxSum() :  GraphMarquee
372        GetMarquee left,bottom
373        if(V_flag == 0)
374                Abort "There is no Marquee"
375        Endif
376        SVAR dispType=root:myGlobals:gDataDisplayType
377        if(cmpstr(dispType,"RealTime")==0)
378                Print "Can't do a BoxSum for a RealTime file"
379                return(1)
380        endif
381        Variable x1,x2,y1,y2
382        x1 = V_left
383        x2 = V_right
384        y1 = V_bottom
385        y2 = V_top
386        KeepSelectionInBounds(x1,x2,y1,y2)
387       
388        String fileStr="",msgStr="Enter a comma-delimited list of run numbers, use dashes for ranges"
389        String type="RAW"
390        Prompt fileStr,msgStr
391        Prompt type,"RAW or Normalized (SAM)",popup,"RAW;SAM;"
392        DoPrompt "Pick the file range",fileStr,type
393        Print "fileStr = ",fileStr
394        printf "(x1,x2) (y1,y2) = (%d,%d) (%d,%d)\r",x1,x2,y1,y2
395       
396        DoBoxSum(fileStr,x1,x2,y1,y2,type)
397       
398        return(0)
399End     
400
401//function that keeps the marquee selection in the range [0,127] inclusive
402// (igor coordinate system)
403// uses pass-by reference!
404//
405// x1 = left
406// x2 = right
407// y1 = bottom
408// y2 = top
409//
410// accepts any detector size
411Function KeepSelectionInBounds(x1,x2,y1,y2)
412        Variable &x1,&x2,&y1,&y2
413       
414        NVAR pixelsX = root:myGlobals:gNPixelsX
415        NVAR pixelsY = root:myGlobals:gNPixelsY
416       
417        //keep selection in-bounds
418        x1 = (round(x1) >= 0) ? round(x1) : 0
419        x2 = (round(x2) <= (pixelsX-1)) ? round(x2) : (pixelsX-1)
420        y1 = (round(y1) >= 0) ? round(y1) : 0
421        y2 = (round(y2) <= (pixelsY-1)) ? round(y2) : (pixelsY-1)
422        return(0)
423End
424
425//testing function, not used
426Function testKeepInBounds(x1,x2,y1,y2)
427        Variable x1,x2,y1,y2
428       
429        KeepSelectionInBounds(x1,x2,y1,y2)
430        Print x1,x2,y1,y2
431        return(0)
432End
433
434// generates a histogram of the data as defined by the marquee. The longer dimension of the marquee
435// becomes the x-axis of the histogram (this may need to be changed for some odd case). Pixel range specified
436// by the marquee is inclusive, and is automatically kept in-bounds
437//
438// The counts over the (short) dimension are averaged, and plotted vs. the pixel position.
439// Pixel position is reported as Detector coordinates (1,128). Counts are whatever the current display
440// happens to be.
441//
442Function SANS_Histogram() :  GraphMarquee
443        GetMarquee left,bottom
444        if(V_flag == 0)
445                Abort "There is no Marquee"
446        endif
447        //
448        Variable count,x1,x2,y1,y2,xwidth,ywidth,vsX=1,xx,yy
449        x1 = V_left
450        x2 = V_right
451        y1 = V_bottom
452        y2 = V_top
453        KeepSelectionInBounds(x1,x2,y1,y2)
454        Print "x1,x2,y1,y2 (det) =",x1+1,x2+1,y1+1,y2+1
455        //determine whether to do x vs y or y vs x
456        xwidth=x2-x1
457        ywidth=y2-y1
458        if(xwidth < ywidth)
459                vsX=0           //sum and graph vs Y
460        endif
461        SVAR cur_folder=root:myGlobals:gDataDisplayType
462        WAVE data=$("root:Packages:NIST:"+cur_folder+":data")           //don't care if it's log or linear scale
463        Make/O/N=(max(xwidth,ywidth)+1) Position,AvgCounts
464        AvgCounts=0
465        //set position wave
466        if(vsX)
467                position=p+x1
468        else
469                position=p+y1
470        endif
471        //convert the position to Detector coordinates
472        position += 1
473       
474        //Compute the histogram (manually)
475        if(vsX)
476                for(xx=x1;xx<=x2;xx+=1)         //outer loop is the "x-axis"
477                        for(yy=y1;yy<=y2;yy+=1)
478                                AvgCounts[xx-x1] += data[xx][yy]
479                        endfor
480                endfor
481                AvgCounts /= (ywidth+1)
482        else
483                for(yy=y1;yy<=y2;yy+=1)
484                        for(xx=x1;xx<=x2;xx+=1)
485                                AvgCounts[yy-y1] += data[xx][yy]
486                        endfor
487                endfor
488                AvgCounts /= (xwidth+1)
489        endif
490        GetMarquee/K            //to keep from drawing the marquee on the new histo graph
491        //draw the graph, or just bring to the front with the new data
492        DoWindow/F SANS_Histo
493        if(V_Flag != 1)
494                Draw_Histo()
495        endif
496       
497        return(0)
498End
499
500//draws the histogram of the 2d data as specified by AvgCounts and Position
501//both wave are assumed to exist in the data folder. The SANS_Histogram() marquee
502//operation is responsible for creating them.
503//
504Function Draw_Histo()
505        Display /W=(197,329,567,461)/K=1 AvgCounts vs Position
506        DoWindow/C SANS_Histo
507        DoWindow/T SANS_Histo,"Histogram"
508        ModifyGraph mode=0,grid=1,mirror=2
509        ModifyGraph rgb=(21845,21845,21845)
510        ModifyGraph standoff=0
511        ModifyGraph hbFill=2
512        ModifyGraph useNegPat=1
513        ModifyGraph usePlusRGB=1
514        ModifyGraph useNegRGB=1
515        ModifyGraph hBarNegFill=2
516        ModifyGraph negRGB=(0,0,65535)
517        SetAxis/A/N=2 left
518        Label left "Counts"
519        Label bottom "Pixel (detector coordinates)"
520End
521
522//function will print marquee coordinates in axis terms, not detector terms
523//since IGOR is [0][127] and detector is (1,128)
524Function PrintMarqueeCoords() :  GraphMarquee
525        GetMarquee left,bottom
526        if(V_flag == 0)
527                Print "There is no Marquee"
528        else
529                Variable count,x1,x2,y1,y2
530                x1 = V_left
531                x2 = V_right
532                y1 = V_bottom
533                y2 = V_top
534                printf "marquee left in bottom axis terms: %g\r",round(V_left)
535                printf "marquee right in bottom axis terms: %g\r",round(V_right)
536                printf "marquee bottom in left axis terms: %g\r",round(V_bottom)
537                printf "marquee top in left axis terms: %g\r",round(V_top)
538                printf "**note that you must add 1 to each axis coordinate to get detector coordinates\r"
539               
540                KeepSelectionInBounds(x1,x2,y1,y2)
541                SVAR cur_folder=root:myGlobals:gDataDisplayType
542                count = SumCountsInBox(x1,x2,y1,y2,cur_folder)
543                Print "counts = "+ num2str(count)
544        endif
545End
Note: See TracBrowser for help on using the repository browser.