Weird wavestats bug..

 
function function1(thresh)
variable thresh

wave wave1,wave2,wave3,wave4,wave5 //preexisting waves in the current experiment

do some stuff
function2(wave1,wave3) // edits wave1 based on information in wave3
do some more stuff
give me nice graphs thanks Igor
end

function2(wave1,wave3)
wave wave1,wave3
wavestats wave1
variable num = V_npntns
wavestats wave3
variable num2 = V_npnts// <-- this is the line that gets highlighted in the debugger, the syntax is identical to the lines above.
do some stuff to wave1 based on num and num2
end

#pragma rtGlobals=3     // Use modern global access method and strict wave access.
function Joint_Hist(thresh)
variable thresh


wave frame, seq_length, intensitych1, intensitych2,ca_conc,amplitudech1,amplitudech2

duplicate/o frame f_key
duplicate/o intensitych1 int_1
duplicate/o intensitych2 int_2
duplicate/o int_2 density
density /= int_1
matrixop/o density = replace(density, NaN, 0)
int_1 = sqrt(intensitych1)
sort int_1,int_1,int_2,f_key,density
join_frames(f_key,seq_length)
wave x_wave = int_1
wave y_wave = f_key

wavestats/q/m=1 y_wave
variable y_dim_max = ceil(V_max)
variable y_dim_min = floor(V_min)
print "y_max is" ,y_dim_max
print "y_min is", y_dim_min
variable y_bins = y_dim_max - y_dim_min + 1
//print y_bins
variable num = V_npnts
wavestats/q/m=1 x_wave
variable x_dim_max = ceil(V_max)
print  "x_max is",x_dim_max
variable x_dim_min = floor(V_min)
print "x_min is", x_dim_min
variable num_bins = 20 // edit this to play with bin resolution
//print num_bins
variable x_delta = (x_dim_max - x_dim_min) / num_bins
print "x_delta is",x_delta
//Build final output waves
make/o/d/n=(num_bins,y_bins) size_hist
make/o/d/n=(num_bins,y_bins,2) density_hist
make/o/d/n=(num_bins,y_bins) density_values
make/o/d/n=(0,3) reject_points
//print "num =",num


//variable thresh = 10  //V_avg + 3*V_sdev
//print "threshold density is",thresh

variable i,j,n
variable x_pos, y_pos
variable max_depth = 0
variable max_errors = 0
for (i=0;i<num;i+=1)
    //print i
    for (j=0;j<num_bins;j+=1)
        if (x_wave[i] - x_dim_min >= (j*x_delta) && x_wave[i] - x_dim_min < (j+1)*x_delta)
            //add a second if-loop to kick out bad data
            if (density[i]<= thresh)
                x_pos = j
                y_pos = y_wave[i] - y_dim_min
                n = size_hist[x_pos][y_pos] // dont iterate n, it corrects the iterate from zero problem
                size_hist[x_pos][y_pos]+= 1
                //check number of binned points and increase matrix depth if needed
                if (n > max_depth)
                    insertpoints/m=2 (max_depth + 1), 1, density_values
                    max_depth +=1
                endif
            density_values[x_pos][y_pos][n] = density[i]       
            else
                max_errors +=1
                insertpoints (max_errors),1, reject_points
                reject_points[max_errors-1][0] = i
                reject_points[max_errors-1][1] = density[i]
                reject_points[max_errors-1][2] = frame[i]
            endif
            break
        endif

    endfor
endfor
// move thru the data matrix and fill in the final histogram
for (i=0;i<num_bins ; i+=1)
    for(j=0;j<y_bins ;j+=1)
        //print "(i , j) =", i,",",j
        n = size_hist[i][j] -1 // if number of counts = 0, the next line duplicates the entire wave.  but the whole wave will be zeros so that is ok.  for now.
        //duplicate/o/r=[i,i][j,j][0,n] density_values temp_wave
        duplicate/o/r=[i,i][j,j][0,max_depth-1] density_values temp_wave
        redimension/n=(max_depth) temp_wave
        wavestats/r=[0,n]/q temp_wave
        density_hist[i][j][0] = V_avg
        density_hist[i][j][1] = V_sdev
    endfor
endfor
newimage/N=density_graph density_hist
SetAxis/W=density_graph/A left
ModifyImage/w=density_graph density_hist ctab= {*,*,Grays,1}
newimage/N=size_graph size_hist
setaxis/W=size_graph/a left
modifyimage/w=size_graph size_hist ctab= {*,*,Grays,1}



end

#pragma rtGlobals=3     // Use modern global access method and strict wave access.
function Join_Frames(frames,seq_lengths)
wave frames, seq_lengths

wavestats/q/m=1 frames
variable i, j
variable num = V_npnts
wavestats/q/m=1 seq_lengths
variable total_seq = V_npnts  // <-- this line gets pointed out in debugger
variable a, b

a=0
for (i=0;i<total_seq;i+=1)
    b= sum(seq_lengths,0,i) - 1
    for (j=0;j<num;j+=1)
        if (frames[j] >= a && frames[j] <=b)
            frames[j] = i
        endif
    endfor
    a = b + 1
endfor
end